!! Parcas_40.f90 is a part of the PACIAE event generator.
!! Copyright (C) 2025 PACIAE Group.
!! PACIAE is licensed under the GNU GPL v2 or later, see LICENSE for details.
!! Open source: https://github.com/ArcsaberHep/PACIAE4
!! Author: Ben-Hao Sa, November 2002 - February 2025.

!> This is the program to deal with the parton cascade (partonic rescattering).

!!                                             By Ben-Hao at CIAE on 19/11/2002
!!                                  Last updated by An-Ke at UiO  on 04/02/2025


        subroutine parcas( time_par )
!!      Deals with the parton cascade (partonic rescattering).
!
!       It was written by Ben-Hao Sa on 19/11/2002.
!       Its input messages are in "PYJETS".
!       Its working block is "PYJETS".
!       Its output messages are in "PYJETS".
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa1/kjp21,non1,bp,iii,neve,nout,nosc
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/sa25/i_inel_proc,i_time_shower,i_diquark,ipad25,para1_1,para1_2
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        common/scatt/ ic, jc, n0
        common/work7/reac(100),crose(100)
!       For the simulation control.
        COMMON/SA1_PY8/ i_mode, i_tune, KF_woDecay(100), &
               KF_proj, KF_targ, win, energy_B, psno, b_min, b_max
        dimension pi00(4), pj00(4)
!       icol: current total number of collision pairs in the collision time list
!       n0: the N before current collision
!       time_par: last collision time
!       icp: the icp-th collision pair
!       tcp: the collision time of the icp-th collision pair
!       i_type: event type. =1, elastic scattering; =2, inelastic; =3, shower.
!       pi00, pj00: four-momenta before collisions


!       Initialization.
        time_par = 0D0
        if( ABS( adj1(1) ) <= 1D-15 ) return
        if( N < 2 ) return
        call reset_eve

!       Creates the (initial) parton-parton collision time list.
        call ctlcre_par
        if( icol == 0 ) return

        !Lei_debug
        ! write(22,*)
        ! write(22,*) "Af call ctlcre_par, iii N, icol =", iii, N, icol
        ! write(22,*) "lc(1-6). tc(1-2), sig, pp2"
        ! do i=1,icol,1
        !     write(22,*) i, ( lc(j,i), j=1,6 ), ( tc(j,i), j=1,2 ) , &
        !                    ( sig_tsmp_Cthetas_p(j,i), j=1,4 )
        ! end do
        ! stop
        ! return
        !Lei_debug


!-------------------------------------------------------------------------------
!-----------------------------   Parton Cascade   ------------------------------
!       Loops over sub-events within an event.
        n_loop = 0
        do while(.true.)
            n0 = N

!       Finds an event (icp) with the minimum time (tcp).
            call find_par( icp, tcp )
            time_par = tcp
            ic     = lc(1,icp)
            jc     = lc(2,icp)
            i_type = lc(6,icp)
            do i=1,4,1
                pi00(i) = P(ic,i)
                pj00(i) = P(jc,i)
            end do

!       Handles a collision (elastic or inelastic) event.
            if( i_type == 1 .OR. i_type == 2 )then
                call collis( icp )

!       Handles a final-state (medium-induced) shower event.
            else
                ! call perform_shower( icp )
            end if

!       Performs classical Newton motion in Lab. system.
            call his_p( ic, jc, pi00, pj00, tcp, i_type )

!       Updates the collision list.
            if( i_time_shower == 0 ) call update_ctl( icp )

            n_loop = n_loop + 1
            if( icol == 0 ) exit
        end do
!-----------------------------   Parton Cascade   ------------------------------
!-------------------------------------------------------------------------------


!       Cross sections statistics.
        do i=1,100,1
            reaci = reac(i)
            if( reaci > 0D0 ) crose(i) = crose(i) / reaci
        end do


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine reset_eve
!!      Initializes the collision time list, etc.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        LOGICAL IS_EXIST, INTO_PRS
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/syspar_p/rsig1,pio,tcut
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        common/scatt/ ic, jc, n0
        common/work7/reac(100),crose(100)
!       For the simulation control.
        COMMON/SA1_PY8/ i_mode, i_tune, KF_woDecay(100), &
               KF_proj, KF_targ, win, energy_B, psno, b_min, b_max
!       ic, jc: current line numbers in "PYJETS" of the collision pair
!       icol: current total number of collision pairs in the collision time list
!       reac and crose: arrays to account for the number and
!        the value of cross sections for corresponding partonic processes
!       taup(i) : formation time of particle i
!       ishp(i)=1, i-th particle will participate in the cascade
!              =0, not


        pio = 3.141592653589793D0
        ic  = 0
        jc  = 0
        n0  = N
        icol  = 0
        reac  = 0D0
        crose = 0D0
        do i=1,N,1
            taup(i) = V(i,4)
            ! Excludes non-partons, historical entries and the
            !  partons that are outside of the simulation time.
            ishp(i) = 0
            KS = K(i,1)
            KF = K(i,2)
            if( .NOT.IS_EXIST(KS,i_mode) .OR. .NOT.INTO_PRS(KF) &
                .OR. V(i,4) > adj1(28) ) cycle
            ishp(i) = 1
        end do


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        LOGICAL FUNCTION INTO_PRS( KF )
!!      Determines whether a particle participates in the partonic rescattering.
        IMPLICIT NONE
        LOGICAL IS_PARTON, IS_DIQUARK
        INTEGER, INTENT(IN) :: KF


        INTO_PRS = IS_PARTON( KF ) .AND. .NOT.IS_DIQUARK( KF )


        RETURN
        END



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine ctlcre_par
!!      Creates the initial collision list.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        integer :: l_coll(6)
        real(kind=8) :: t_coll, t_coll_min
        real(kind=8) :: sig_t_theta_p(4)
        logical :: succeed


        ! Time resolution.
        dddt = adj1(19)
        ! The collision pair counter.
        icol = 1
        do i = 1, N-1, 1
            if( ishp(i) == 0 ) cycle
            ! Successful flag.
            succeed = .false.
            ! A number large enough to determine the minumum time.
            t_coll_min = 1D30
            ixc = 0
            jxc = 0
            loop_j: do j = i+1, N, 1
                if( ishp(j) == 0 ) cycle
                call coij_p( i, j, i_fail, l_coll, t_coll, sig_t_theta_p )
                if( i_fail == 1 ) cycle
                ! Imposes the time resolution constraint.
                if( t_coll < 1D-10 ) cycle
                do j1 = 1, icol-1, 1
                    if( ABS( tc(1,j1) - t_coll ) < dddt ) cycle loop_j
                end do
                ! Chooses the smallest time for 'i'-cycle.
                if( t_coll < t_coll_min )then
                    lc(1,icol) = l_coll(1)
                    lc(2,icol) = l_coll(2)
                    lc(3,icol) = l_coll(3)
                    lc(4,icol) = l_coll(4)
                    lc(5,icol) = l_coll(5)
                    lc(6,icol) = l_coll(6)
                    tc(1,icol) = t_coll
                    tc(2,icol) = t_coll
                    sig_tsmp_Cthetas_p(1,icol) = sig_t_theta_p(1)
                    sig_tsmp_Cthetas_p(2,icol) = sig_t_theta_p(2)
                    sig_tsmp_Cthetas_p(3,icol) = sig_t_theta_p(3)
                    sig_tsmp_Cthetas_p(4,icol) = sig_t_theta_p(4)
                    t_coll_min = t_coll
                    ixc = l_coll(1)
                    jxc = l_coll(2)
                end if
                succeed = .true.
            end do loop_j
            ! Keeps the one with the smallest time from pairs including i or j.
            if( succeed )then
                n_jump_out = 0
                j1 = 1
                do while(.true.)
                    if( j1 > icol-1 ) exit
                    iic = lc(1,j1)
                    jjc = lc(2,j1)
                    if(       ixc /= iic .AND. ixc /= jjc &
                        .AND. jxc /= iic .AND. jxc /= jjc )then
                        j1  = j1 + 1
                        cycle
                    end if
                    ttc = tc(1,j1)
                    ! Thows away the pair with larger time.
                    if( ttc > t_coll_min )then
                        k_begin = j1
                        n_jump_out = n_jump_out + 1
                    else
                        k_begin = icol
                        n_jump_out = 2
                    end if
                    icol = icol - 1
                    do k1 = k_begin, icol, 1
                        lc( 1, k1 ) = lc( 1, k1+1 )
                        lc( 2, k1 ) = lc( 2, k1+1 )
                        lc( 3, k1 ) = lc( 3, k1+1 )
                        lc( 4, k1 ) = lc( 4, k1+1 )
                        lc( 5, k1 ) = lc( 5, k1+1 )
                        lc( 6, k1 ) = lc( 6, k1+1 )
                        tc( 1, k1 ) = tc( 1, k1+1 )
                        tc( 2, k1 ) = tc( 2, k1+1 )
                        sig_tsmp_Cthetas_p(1,k1) = sig_tsmp_Cthetas_p( 1, k1+1 )
                        sig_tsmp_Cthetas_p(2,k1) = sig_tsmp_Cthetas_p( 2, k1+1 )
                        sig_tsmp_Cthetas_p(3,k1) = sig_tsmp_Cthetas_p( 3, k1+1 )
                        sig_tsmp_Cthetas_p(4,k1) = sig_tsmp_Cthetas_p( 4, k1+1 )
                    end do
                    if( n_jump_out == 2 ) exit
                end do
                ! For the next time finding.
                icol = icol + 1
            end if
        end do
        ! Deducts 1 because the counter begins with 1 not 0.
        icol = icol - 1


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine find_par( icp, tcp )
!!      Finds out the binary collision (icp) with minimum colli. time (tcp).
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (MCLIS=280000)
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol


        icp = MINLOC( tc( 1, 1:icol ), DIM=1 )
        tcp = MINVAL( tc( 1, 1:icol ), DIM=1 )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine update_ctl( icp )
!!      Throws away old collission pairs and updates the collision time list.
!       "update_ctl" = "update_ctlm" + "update_ctln"
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/scatt/ ic, jc, n0
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        integer, intent(in) :: icp
        integer :: l_coll(6)
        real(kind=8) :: t_coll, t_coll_min
        real(kind=8) :: sig_t_theta_p(4)
        logical :: succeed


!       Throws away the currently scattered pair.
        icol = icol - 1
        do k1 = icp, icol, 1
            lc( 1, k1 ) = lc( 1, k1+1 )
            lc( 2, k1 ) = lc( 2, k1+1 )
            lc( 3, k1 ) = lc( 3, k1+1 )
            lc( 4, k1 ) = lc( 4, k1+1 )
            lc( 5, k1 ) = lc( 5, k1+1 )
            lc( 6, k1 ) = lc( 6, k1+1 )
            tc( 1, k1 ) = tc( 1, k1+1 )
            tc( 2, k1 ) = tc( 2, k1+1 )
            sig_tsmp_Cthetas_p(1,k1) = sig_tsmp_Cthetas_p( 1, k1+1 )
            sig_tsmp_Cthetas_p(2,k1) = sig_tsmp_Cthetas_p( 2, k1+1 )
            sig_tsmp_Cthetas_p(3,k1) = sig_tsmp_Cthetas_p( 3, k1+1 )
            sig_tsmp_Cthetas_p(4,k1) = sig_tsmp_Cthetas_p( 4, k1+1 )
        end do

!       Loops over ic,jc (new) and old partons (i.e. constructs colli. pair
!        by partons, one of which is ic/jc and another one is in parton list).

        ! Time resolution.
        dddt = adj1(19)
        ! The collision pair counter (+1 avoids 0 icol).
        icol = icol + 1
        do ii=1,2,1
            if( ii == 1 ) i = ic
            if( ii == 2 ) i = jc
            if( ishp(i) == 0 ) cycle
            ! Successful flag.
            succeed = .false.
            ! A number large enough to determine the minumum time.
            t_coll_min = 1D30
            ixc = 0
            jxc = 0
            loop_j: do j = 1, n0, 1
                if( ishp(j) == 0 ) cycle
                ! Avoids two particles colliding immediately after a collision.
                if( j == ic .OR. j == jc ) cycle
                call coij_p( i, j, i_fail, l_coll, t_coll, sig_t_theta_p )
                if( i_fail == 1 ) cycle
                ! Imposes the time resolution constraint.
                if( t_coll < 1D-10 ) cycle
                do j1 = 1, icol-1, 1
                    if( ABS( tc(1,j1) - t_coll ) < dddt ) cycle loop_j
                end do
                ! Chooses the smallest time for 'i'-cycle.
                if( t_coll < t_coll_min )then
                    lc(1,icol) = l_coll(1)
                    lc(2,icol) = l_coll(2)
                    lc(3,icol) = l_coll(3)
                    lc(4,icol) = l_coll(4)
                    lc(5,icol) = l_coll(5)
                    lc(6,icol) = l_coll(6)
                    tc(1,icol) = t_coll
                    tc(2,icol) = t_coll
                    sig_tsmp_Cthetas_p(1,icol) = sig_t_theta_p(1)
                    sig_tsmp_Cthetas_p(2,icol) = sig_t_theta_p(2)
                    sig_tsmp_Cthetas_p(3,icol) = sig_t_theta_p(3)
                    sig_tsmp_Cthetas_p(4,icol) = sig_t_theta_p(4)
                    t_coll_min = t_coll
                    ixc = l_coll(1)
                    jxc = l_coll(2)
                end if
                succeed = .true.
            end do loop_j
            ! Keeps the one with the smallest time from pairs including i or j.
            if( succeed )then
                n_jump_out = 0
                j1 = 1
                do while(.true.)
                    if( j1 > icol-1 ) exit
                    iic = lc(1,j1)
                    jjc = lc(2,j1)
                    if(       ixc /= iic .AND. ixc /= jjc &
                        .AND. jxc /= iic .AND. jxc /= jjc )then
                        j1  = j1 + 1
                        cycle
                    end if
                    ttc = tc(1,j1)
                    ! Thows away the pair with larger time.
                    if( ttc > t_coll_min )then
                        k_begin = j1
                        n_jump_out = 2
                    else
                        k_begin = icol
                        n_jump_out = 2
                    end if
                    icol = icol - 1
                    do k1 = k_begin, icol, 1
                        lc( 1, k1 ) = lc( 1, k1+1 )
                        lc( 2, k1 ) = lc( 2, k1+1 )
                        lc( 3, k1 ) = lc( 3, k1+1 )
                        lc( 4, k1 ) = lc( 4, k1+1 )
                        lc( 5, k1 ) = lc( 5, k1+1 )
                        lc( 6, k1 ) = lc( 6, k1+1 )
                        tc( 1, k1 ) = tc( 1, k1+1 )
                        tc( 2, k1 ) = tc( 2, k1+1 )
                        sig_tsmp_Cthetas_p(1,k1) = sig_tsmp_Cthetas_p( 1, k1+1 )
                        sig_tsmp_Cthetas_p(2,k1) = sig_tsmp_Cthetas_p( 2, k1+1 )
                        sig_tsmp_Cthetas_p(3,k1) = sig_tsmp_Cthetas_p( 3, k1+1 )
                        sig_tsmp_Cthetas_p(4,k1) = sig_tsmp_Cthetas_p( 4, k1+1 )
                    end do
                    if( n_jump_out == 2 ) exit
                end do
                ! For the next time finding.
                icol = icol + 1
            end if
        end do
        ! Deducts 1 because the counter begins with +1 not 0.
        icol = icol - 1
        n0 = N


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine update_ctlm( icp )
!!      Throws away old collission pairs.
!       "update_ctl" = "update_ctlm" + "update_ctln"
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (MCLIS=280000)
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        integer, intent(in) :: icp


!       Throws away the currently scattered pair.
        icol = icol - 1
        do k1 = icp, icol, 1
            lc( 1, k1 ) = lc( 1, k1+1 )
            lc( 2, k1 ) = lc( 2, k1+1 )
            lc( 3, k1 ) = lc( 3, k1+1 )
            lc( 4, k1 ) = lc( 4, k1+1 )
            lc( 5, k1 ) = lc( 5, k1+1 )
            lc( 6, k1 ) = lc( 6, k1+1 )
            tc( 1, k1 ) = tc( 1, k1+1 )
            tc( 2, k1 ) = tc( 2, k1+1 )
            sig_tsmp_Cthetas_p(1,k1) = sig_tsmp_Cthetas_p( 1, k1+1 )
            sig_tsmp_Cthetas_p(2,k1) = sig_tsmp_Cthetas_p( 2, k1+1 )
            sig_tsmp_Cthetas_p(3,k1) = sig_tsmp_Cthetas_p( 3, k1+1 )
            sig_tsmp_Cthetas_p(4,k1) = sig_tsmp_Cthetas_p( 4, k1+1 )
        end do


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine update_ctln
!!      Updates the collision time list.
!       "update_ctl" = "update_ctlm" + "update_ctln"
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/scatt/ ic, jc, n0
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        integer :: l_coll(6)
        real(kind=8) :: t_coll, t_coll_min
        real(kind=8) :: sig_t_theta_p(4)
        logical :: succeed


!       Loops over ic,jc (new) and old partons (i.e. constructs colli. pair
!        by partons, one of which is ic/jc and another one is in parton list).

        ! Time resolution.
        dddt = adj1(19)
        ! The collision pair counter (+1 avoids 0 icol).
        icol = icol + 1
        do ii=1,2,1
            if( ii == 1 ) i = ic
            if( ii == 2 ) i = jc
            if( ishp(i) == 0 ) cycle
            ! Successful flag.
            succeed = .false.
            ! A number large enough to determine the minumum time.
            t_coll_min = 1D30
            ixc = 0
            jxc = 0
            loop_j: do j = 1, n0, 1
                if( ishp(j) == 0 ) cycle
                ! Avoids two particles colliding immediately after a collision.
                if( j == ic .OR. j == jc ) cycle
                call coij_p( i, j, i_fail, l_coll, t_coll, sig_t_theta_p )
                if( i_fail == 1 ) cycle
                ! Imposes the time resolution constraint.
                if( t_coll < 1D-10 ) cycle
                do j1 = 1, icol-1, 1
                    if( ABS( tc(1,j1) - t_coll ) < dddt ) cycle loop_j
                end do
                ! Chooses the smallest time for 'i'-cycle.
                if( t_coll < t_coll_min )then
                    lc(1,icol) = l_coll(1)
                    lc(2,icol) = l_coll(2)
                    lc(3,icol) = l_coll(3)
                    lc(4,icol) = l_coll(4)
                    lc(5,icol) = l_coll(5)
                    lc(6,icol) = l_coll(6)
                    tc(1,icol) = t_coll
                    tc(2,icol) = t_coll
                    sig_tsmp_Cthetas_p(1,icol) = sig_t_theta_p(1)
                    sig_tsmp_Cthetas_p(2,icol) = sig_t_theta_p(2)
                    sig_tsmp_Cthetas_p(3,icol) = sig_t_theta_p(3)
                    sig_tsmp_Cthetas_p(4,icol) = sig_t_theta_p(4)
                    t_coll_min = t_coll
                    ixc = l_coll(1)
                    jxc = l_coll(2)
                end if
                succeed = .true.
            end do loop_j
            ! Keeps the one with the smallest time from pairs including i or j.
            if( succeed )then
                n_jump_out = 0
                j1 = 1
                do while(.true.)
                    if( j1 > icol-1 ) exit
                    iic = lc(1,j1)
                    jjc = lc(2,j1)
                    if(       ixc /= iic .AND. ixc /= jjc &
                        .AND. jxc /= iic .AND. jxc /= jjc )then
                        j1  = j1 + 1
                        cycle
                    end if
                    ttc = tc(1,j1)
                    ! Thows away the pair with larger time.
                    if( ttc > t_coll_min )then
                        k_begin = j1
                        n_jump_out = 2
                    else
                        k_begin = icol
                        n_jump_out = 2
                    end if
                    icol = icol - 1
                    do k1 = k_begin, icol, 1
                        lc( 1, k1 ) = lc( 1, k1+1 )
                        lc( 2, k1 ) = lc( 2, k1+1 )
                        lc( 3, k1 ) = lc( 3, k1+1 )
                        lc( 4, k1 ) = lc( 4, k1+1 )
                        lc( 5, k1 ) = lc( 5, k1+1 )
                        lc( 6, k1 ) = lc( 6, k1+1 )
                        tc( 1, k1 ) = tc( 1, k1+1 )
                        tc( 2, k1 ) = tc( 2, k1+1 )
                        sig_tsmp_Cthetas_p(1,k1) = sig_tsmp_Cthetas_p( 1, k1+1 )
                        sig_tsmp_Cthetas_p(2,k1) = sig_tsmp_Cthetas_p( 2, k1+1 )
                        sig_tsmp_Cthetas_p(3,k1) = sig_tsmp_Cthetas_p( 3, k1+1 )
                        sig_tsmp_Cthetas_p(4,k1) = sig_tsmp_Cthetas_p( 4, k1+1 )
                    end do
                    if( n_jump_out == 2 ) exit
                end do
                ! For the next time finding.
                icol = icol + 1
            end if
        end do
        ! Deducts 1 because the counter begins with +1 not 0.
        icol = icol - 1
        n0 = N


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine coij_p( ic, jc, i_fail, l_coll, t_coll, sig_t_theta_p )
!!      Calculates the collision time and kinematics of two particles.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/syspar_p/rsig1,pio,tcut
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        integer, intent(in) :: ic, jc
        integer :: l_coll(6)
        real(kind=8) :: t_coll
        real(kind=8) :: sig_t_theta_p(4)
        dimension ri(4), rj(4), ric(4), rjc(4)
        dimension pi(4), pj(4), pic(4), pjc(4), pij(4)
        dimension dr(3), db(3), vi(3), vj(3)
        dimension b(3)


        i_fail = 1
        l_coll = 0
        t_coll = 0D0
        sig_t_theta_p = 0D0
        t_max  = adj1(28)

        pi(4)  = P(ic,4)
        pj(4)  = P(jc,4)
        do k1=1,3,1
            pi(k1) = P(ic,k1)
            pj(k1) = P(jc,k1)
            b(k1)  = ( pi(k1) + pj(k1) ) / ( pi(4) + pj(4) )
        end do
        pic = pi
        pjc = pj
        pij = pi + pj
!       Invariant mass squared.
        eiej2 = pij(4)**2 - pij(1)**2 - pij(2)**2 - pij(3)**2
!       Inserts the energy cut.
        dm1  = P(ic,5)
        dm2  = P(jc,5)
        ecut = SQRT(eiej2) - dm1 - dm2
        if( ecut <= 0D0 ) return
        do n1=1,4,1
            ri(n1)  = V(ic,n1)
            rj(n1)  = V(jc,n1)
        end do
        ric = ri
        rjc = rj
        KF1 = K(ic,2)
        KF2 = K(jc,2)
        call fsig( KF1, KF2, KF3, KF4, dm1, dm2, dm3, dm4, eiej2, &
                   sig, tsmp, cos_theta_s, i_proc )
        if( sig <= 0D0 ) return
        rsig1 = SQRT( sig / pio )

!       Boosts into the CM frame.
        ! Momentum.
        call lorntz( 0, b, pic, pjc )
        ! Position.
        call lorntz( 0, b, ric, rjc )

!       Finds the minimum distance and the time.
        rb = 0D0
        bb = 0D0
        rtai = 0D0
        do k1=1,3,1
            vi(k1) = pic(k1) / pic(4)
            vj(k1) = pjc(k1) / pjc(4)
        end do
        do k1=1,3,1
            dr(k1) = ric(k1) - rjc(k1) - ( vi(k1)*ric(4) - vj(k1)*rjc(4) )
            db(k1) = vi(k1)  - vj(k1)
            rb = rb + dr(k1)*db(k1)
            bb = db(k1)**2 + bb
        end do
        if( bb <= 1D-10 ) return
        tcol = 0D0 - rb / bb
!       The collision should happen in the future of the CM frame.
        if( tcol < MAX( ric(4), rjc(4) ) ) return
        do ik=1,3,1
            dr(ik) = dr(ik) + tcol * db(ik)
            rtai = rtai + dr(ik) * dr(ik)
        end do
        sg = rtai
        dmin = SQRT(sg)
        if( dmin > rsig1 ) return
!       Moves along the Newton trajectory in CMS.
        do ik=1,3,1
            ric(ik) = ric(ik) + vi(ik) * ( tcol - ric(4) )
            rjc(ik) = rjc(ik) + vj(ik) * ( tcol - rjc(4) )
        end do
        ric(4) = tcol
        rjc(4) = tcol

!       Transforms back to the Lab frame (causality violation).
        call lorntz( 1, b, ric, rjc )
!       Chooses the min one for the causality violation.
        t_coll = MIN( ric(4), rjc(4) )
!       The collision should happen in the future of the Lab frame.
        if( t_coll < MAX( ri(4), rj(4), 0D0 ) ) return
!       Max time constraint (default 10000 fm/c).
        if( t_coll > t_max ) return

!       Kinematics of inelastic scatterings.
        i_proc1 = MOD( i_proc, 10 )
        if( i_proc1 == 2 .OR. i_proc1 == 6 .OR. i_proc1 ==7 )then
            i_type = 2
            pp2 = ( eiej2 - (dm3 + dm4)**2 ) * ( eiej2 - (dm3 - dm4)**2 ) &
                / 4D0 / eiej2
!       Kinematics of elastic scatterings.
        else
            i_type = 1
            pp2 = pic(1)**2 + pic(2)**2 + pic(3)**2
        end if
!       Too small/zero momentum.
        if( pp2 < 1D-20 ) return

!       Succeeded.
        l_coll(1) = ic
        l_coll(2) = jc
        l_coll(3) = KF3
        l_coll(4) = KF4
        l_coll(5) = i_proc
        l_coll(6) = i_type
        sig_t_theta_p(1) = sig
        sig_t_theta_p(2) = tsmp
        sig_t_theta_p(3) = cos_theta_s
        sig_t_theta_p(4) = SQRT( pp2 )

        i_fail = 0


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_dot(v1,v2)
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        dimension v1(4),v2(4)


        func_dot = v1(4)*v2(4) - v1(1)*v2(1) - v1(2)*v2(2) - v1(3)*v2(3)


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine fsig( kf1, kf2, kf3, kf4, m1, m2, m3, m4, s, &
                         sig, tsmp, cos_theta_s, i_proc  )
!!      Calculates the total cross section, decides the type of reaction,
!!       samples the t value and the scattering angle cos(theta_s).
!
!       The masses of the heavy c and b are retained (massive) in any case, and
!        other masses of light quarks might treated to zero (zero mass
!        approximation in kinematics treatment, massless).
!
!       kf1, kf2: kf codes of the incoming pair.
!       kf3, kf4: kf codes of the outgoing pair.
!       m1, m2: masses of the incoming pair.
!       m3, m4: masses of the outgoing pair.
!       s: the invariant mass squared of the colliding pair.
!       sig: the total cross section of parton kf1 colliding with kf2.
!       tsmp: the momentum transfer value sampled.
!       cos_theta_s: the scattering angle COS(theta_s).
!       i_proc: the internal processes number.
!
!       i_proc:
!        In the following internal order, we do not distinguish between q and Q.
!           1: g     + g     -> g     + g
!           2: g     + g     -> q1    + q1bar   (inelastic)
!           3: q     + g     -> q     + g
!              qbar  + g     -> qbar  + g
!           4: q1    + q1    -> q1    + q1
!              q1bar + q1bar -> q1bar + q1bar
!           5: q1    + q1bar -> q1    + q1bar
!           6: q1    + q1bar -> q2    + q2bar   (inelastic)
!           7: q1    + q1bar -> g     + g       (inelastic)
!           8: q1    + q2    -> q1    + q2
!              q1    + q2bar -> q1    + q2bar
!              q1bar + q2bar -> q1bar + q2bar
!#TODO(Lei20241113): Heavy onia gluon dissociation.
!          xx: hOnia + g  -> hQ + hQbar       (dissociation)
!#TODO(Lei20241113): Heavy onia gluon regeneration.
!          xx: hQ + hQbar -> hOnia + g        (regeneration)
!         +10: with one heavy quark
!         +20: with two heavy quark
!
!       References:
!        Massless:
!            B.L. Combridge et al., Phys. Lett. B 70 (1977) 234;
!            Bin Zhang, Comput.Phys.Commun. 109 (1998) 193-206;
!            Jussi Auvinen et al., Phys.Rev.C 82 (2010) 024906;
!            Ben-Hao Sa, Comput.Phys.Commun. 183 (2012) 333–346.
!        Massive:
!            B.L. Combridge Nucl. Phys. B 151 (1979) 429-456;
!            Hamza Berrehrah et al., Phys.Rev.C 89 (2014) 5, 054901.
!            V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
!        t sampling:
!            Klaus Geiger et al., Nucl.Phys.B 369 (1992) 600-654;
!            Ben-Hao Sa et al., Comput.Phys.Commun. 183 (2012) 333–346.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/sa25/i_inel_proc,i_time_shower,i_diquark,ipad25,para1_1,para1_2
        common/sa33/smadel,ecce,secce,parecc,iparres
        integer, intent(in) :: kf1, kf2
        real(kind=8), intent(in) :: m1, m2, s
        real(kind=8) :: m3, m4
!       Note the following logical type.
        logical :: l_gg
        logical :: l_qg, l_qbarg
        logical :: l_q1q1, l_q1barq1bar
        logical :: l_q1q1bar
        logical :: l_q1q2, l_q1q2bar, l_q1barq2bar
        logical :: l_ghOnia
        real(kind=8) :: sigma_relative( 0:7 )
        real(kind=8) :: t_0_relative( 0:7 ),   t_pi_relative( 0:7 ), &
                        t_max_relative( 0:7 ), t_min_relative( 0:7 )
        real(kind=8) :: alpha_s, Lambda_QCD2, K_factor
        integer :: KF_outgoing( 2, 0:7 )


        sig  = 0D0
        tsmp = 0D0
        cos_theta_s = 0D0
        i_proc = 0
        kf3  = 0
        kf4  = 0
        dm3  = 0D0
        dm4  = 0D0
        m3   = 0D0
        m4   = 0D0
        if( s <= 0D0 ) return


!-------------------------------------------------------------------------------
!---------------------------   Process Classifying   ---------------------------
!       Classifies the incident partons.
!       kf: 1 ~ 5 -- d, u, s, c, b; 21 -- g
!           443 -- J/psi, 100443 -- Psi';
!           553 -- Upsilon, 100553 -- Upsilon';
        ikf1 = ABS(kf1)
        ikf2 = ABS(kf2)

        !Lei_debug
        if( ikf1 >= 4 .AND. ikf1 <= 5 .AND. ikf2 >= 4 .AND. ikf2 <= 5 ) return
        !Lei_debug

        l_gg = .false.
        l_qg = .false.
        l_qbarg = .false.
        l_q1q1  = .false.
        l_q1barq1bar = .false.
        l_q1q1bar    = .false.
        l_q1q2       = .false.
        l_q1q2bar    = .false.
        l_q1barq2bar = .false.
        l_ghOnia     = .false.
        ! g + g ->
        if( kf1 == 21 .AND. kf2 == 21 )then
            l_gg = .true.
        ! X + g -> X + g / X
        else if( kf1 == 21 .OR. kf2 == 21 )then
            ! q + g -> q + g
            if(      ( kf1 >= 1 .AND. kf1 <= 8 ) &
                .OR. ( kf2 >= 1 .AND. kf2 <= 8 ) )then
                    l_qg = .true.
            ! qbar + g -> qbar + g
            else if( ( kf1 <= -1 .AND. kf1 >= -8 ) &
                .OR. ( kf2 <= -1 .AND. kf2 >= -8 ) )then
                    l_qbarg = .true.
            ! g + hOnia -> hQ + hQbar
            ! Only J/psi, Psi', Upsilon and and Upsilon'.
            else if( ( ikf1 == 443 .OR. ikf1 == 100443   &
                .OR.   ikf1 == 553 .OR. ikf1 == 100553 ) &
                .OR. ( ikf2 == 443 .OR. ikf2 == 100443   &
                .OR.   ikf2 == 553 .OR. ikf2 == 100553 ) )then
                    l_ghOnia = .true.
            end if
        ! q + q ->
        else if( ikf1 >= 1 .AND. ikf1 <= 8 .AND. ikf2 >= 1 .AND. ikf2 <= 8 )then
            ! q1 + q1 -> q1 + q1
            if( kf1 == kf2 .AND. kf1 > 0 )then
                l_q1q1 = .true.
            ! q1bar + q1bar -> q1bar + q1bar
            else if( kf1 == kf2 .AND. kf1 < 0 )then
                l_q1barq1bar = .true.
            ! q1 + q1bar ->
            else if( ikf1 == ikf2 .AND. kf1*kf2 < 0 )then
                l_q1q1bar = .true.
            ! q1 + q2 -> q1 + q2
            else if( ikf1 /= ikf2 .AND. kf1 > 0 .AND. kf2 > 0 )then
                l_q1q2 = .true.
            ! q1 + q2bar -> q1 + q2bar
            else if( ikf1 /= ikf2 .AND. kf1*kf2 < 0 )then
                l_q1q2bar = .true.
            ! q1bar + q2bar -> q1bar + q2bar
            else if( ikf1 /= ikf2 .AND. kf1 < 0 .AND. kf2 < 0 )then
                l_q1barq2bar = .true.
            end if
        else
            write(*,*) "Warning! fsig: wrong KFs ! kf1, kf2 =", kf1, kf2
            return
        end if
!---------------------------   Process Classifying   ---------------------------
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!--------------------------   Variable Initializing   --------------------------
!       Relative cross sections among different outgoing partons in the
!        inelastic processes.
        sigma_relative = 0D0
        t_0_relative   = 0D0
        t_pi_relative  = 0D0
        t_max_relative = 0D0
        t_min_relative = 0D0
        KF_outgoing    = 0
!       Switch of inelastic processes.
        i_inelastic = iparres
!       Selection of inelastic processes when i_inelastic = 1.
!       i_inel_proc: = 1, with only light quarks related inelastic processes
!                    = 2, with both light and heavy quarks related inel. proc.
        n_flavor_sample = 0
        if( i_inelastic == 1 )then
            if( i_inel_proc == 1 )then
                n_flavor_sample = 3
            else if( i_inel_proc == 2 )then
                n_flavor_sample = 5
            end if
        end if
!       Selection of the form of differential cross section in the LO-pQCD.
!       i_LOpQCD:
!        =1, small angle and zero-mass approximations for light quarks in
!            elastic parton-parton cross sections with the regulator.
!        =2, full LO-pQCD and massive forms for parton-parton cross sections
!            with the regulator.
!        =3, small angle and zero-mass approximations for light quarks in
!            elastic parton-parton cross sections
!            with cutoffs in the integral limits.
!        =4, full LO-pQCD and massive forms for parton-parton cross sections
!            with cutoffs in the integral limits.
!       The inelastic processes always utilize non-zero-mass full LO-QCD forms.
        i_LOpQCD = INT( adj1(20) )
        dm1 = m1
        dm2 = m2
        ! Zero-mass approximation for light quarks.
        if( i_LOpQCD == 1 .OR. i_LOpQCD == 3 )then
            if( ikf1 < 4 ) dm1 = 0D0
            if( ikf2 < 4 ) dm2 = 0D0
        end if
!       Phenomenological high-order correction factor K.
        K_factor = adj1(1)
!       Invariant mass of the colliding pair.
        sqrt_s = SQRT( s )
        ! Basic threshold energy.
        if( sqrt_s < ( m1 + m2 ) ) return
!       Effective coupling constant.
        ! Constant alpha_s.
        alpha_s  = adj1(2)
        alpha_s0 = adj1(2)
        ! Renormalization scale mu_R = Q.
        Q2 = s - dm1*dm1 - dm2*dm2
        ! Number of active quarks.
        Nf = 3
        if( Q2 > amass(4)*amass(4) ) Nf = 4
        if( Q2 > amass(5)*amass(5) ) Nf = 5
        ! QCD scale.
        Lambda_QCD2 = adj1(25)*adj1(25)
        ! Running alpha_s.
        if( alpha_s0 < 1D-10 ) alpha_s = func_alpha_s( Nf, Q2, Lambda_QCD2 )
        if( alpha_s  < 1D-10 ) return
!       t_cut: the cutoff used to regulate the 't'/'u' divergence.
        ! Constant cutoff (screening masses).
        t_cut  = adj1(3)
        t_cut0 = adj1(3)
        ! Gluon screening mass.
        dmD2   = t_cut
        ! Quark medium mass.
        dmq2   = t_cut
        ! Dynamic screening masses.
        if( t_cut0 < 1D-10 )then
            T_effictive = ABS( adj1(3) )
            ! Boltzmann or Bose-Einstein & Fermi-Dirac statistics in screening.
            ! i_statistic = 0
            i_statistic = 1
            ! Gluon screening mass.
            dmD2  = func_mD2( alpha_s, Nf, T_effictive, i_statistic )
            t_cut = dmD2
            ! Quark medium mass.
            dmq2  = func_mq2( alpha_s, T_effictive, i_statistic )
        end if
!       Momentum transfer limits (scattering angles = 0 & pi) of massless proc.
        t_0  = -s
        t_pi = 0D0
!       Note that the integral variable t is negative.
!--------------------------   Variable Initializing   --------------------------
!-------------------------------------------------------------------------------


!       Calculates the total cross section for each colliding pair.


!-------------------------------------------------------------------------------
!------------------------   Cross Section Calculating   ------------------------
!       g + g ->
        if( l_gg )then
!       1. Elastic process: g + g -> g + g.
            i_proc = 1
            dm3    = 0D0
            dm4    = 0D0
            t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
            t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
            t_pi_relative(0) = t_min
            t_0_relative(0)  = t_max
            if( i_LOpQCD > 2 )then
                t_max = t_max - t_cut
                if( i_LOpQCD == 4 ) t_min = t_min + t_cut
            end if
            t_min_relative(0) = t_min
            t_max_relative(0) = t_max
            sig_out = get_sigma( t_min, t_max, t_cut, s, dm3, dm4, &
                                 alpha_s, K_factor, i_proc, i_LOpQCD )
            sigma_relative(0)   = sig_out
            KF_outgoing( 1, 0 ) = 21
            KF_outgoing( 2, 0 ) = 21
!       2. Inelastic process: g + g -> q + qbar.
            if( i_inelastic == 1 )then
                i_proc = 2
                do i_flavor = 1, n_flavor_sample, 1
                    ! KF = 1 : d; KF = 2 : u.
                    i_flavor_out = i_flavor
                    if( i_flavor == 1 ) i_flavor_out = 2
                    if( i_flavor == 2 ) i_flavor_out = 1
                    ! Threshold energy.
                    if( sqrt_s >= 2D0*amass( i_flavor_out ) )then
                        dm3    = amass( i_flavor_out )
                        dm4    = amass( i_flavor_out )
                        t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
                        t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
                        t_pi_relative( i_flavor )  = t_min
                        t_0_relative(  i_flavor )  = t_max
                        t_max = t_max - t_cut
                        t_min = t_min + t_cut
                        t_min_relative( i_flavor ) = t_min
                        t_max_relative( i_flavor ) = t_max
                        sig_out = get_sigma( t_min, t_max, t_cut, s,      &
                                             dm3, dm4, alpha_s, K_factor, &
                                             i_proc, i_LOpQCD )
                        sigma_relative( i_flavor ) = sig_out
                        KF_outgoing( 1, i_flavor ) =  i_flavor_out
                        KF_outgoing( 2, i_flavor ) = -i_flavor_out
                    end if
                end do
            end if
!           Calculates the total corss section and determines the process to be
!            happened according to the relative cross section.
            sig = 0D0
            do i = 0, n_flavor_sample, 1
                sig = sig + sigma_relative(i)
            end do
            rand_num = PYR(1)
            prob_low = 0D0
            prob_upp = 0D0
            do i = 0, n_flavor_sample, 1
                prob_upp = prob_upp + sigma_relative(i) / sig
                if( rand_num > prob_low .AND. rand_num < prob_upp ) exit
                prob_low = prob_upp
            end do
            i_proc = 1
            if( i > 0 ) i_proc = 2
            kf3 = KF_outgoing( 1, i )
            kf4 = KF_outgoing( 2, i )
            dm3 = amass( kf3 )
            dm4 = amass( kf4 )
            t_min = t_min_relative( i )
            t_max = t_max_relative( i )
            t_pi  = t_pi_relative( i )
            t_0   = t_0_relative( i )
!       3. Elastic process: q + g -> q + g or qbar + g -> qbar + g.
        else if( l_qg .OR. l_qbarg )then
            i_proc = 3
            dm3    = dm1
            dm4    = dm2
            t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
            t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
            t_pi   = t_min
            t_0    = t_max
            if( i_LOpQCD > 2 )then
                t_max = t_max - t_cut
                if( i_LOpQCD == 4 ) t_min = t_min + t_cut
            end if
            sig = get_sigma( t_min, t_max, t_cut, s, dm3, dm4, &
                             alpha_s, K_factor, i_proc, i_LOpQCD )
            kf3 = kf1
            kf4 = kf2
!       4. Elastic process: q1 + q1 -> q1 + q1, q1bar + q1bar -> q1bar + q1bar.
        else if( l_q1q1 .OR. l_q1barq1bar )then
            i_proc = 4
            dm3    = dm1
            dm4    = dm2
            t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
            t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
            t_pi   = t_min
            t_0    = t_max
            if( i_LOpQCD > 2 )then
                t_max = t_max - t_cut
                if( i_LOpQCD == 4 ) t_min = t_min + t_cut
            end if
            sig = get_sigma( t_min, t_max, t_cut, s, dm3, dm4, &
                             alpha_s, K_factor, i_proc, i_LOpQCD )
            kf3 = kf1
            kf4 = kf2
!       q1 + q1bar ->
        else if( l_q1q1bar )then
!       5. Elastic process: q1 + q1bar -> q1 + q1bar.
            i_proc = 5
            dm3    = dm1
            dm4    = dm2
            t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
            t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
            t_pi_relative(0) = t_min
            t_0_relative(0)  = t_max
            if( i_LOpQCD > 2 )then
                t_max = t_max - t_cut
            end if
            t_min_relative(0) = t_min
            t_max_relative(0) = t_max
            sig_out = get_sigma( t_min, t_max, t_cut, s, dm3, dm4, &
                                 alpha_s, K_factor, i_proc, i_LOpQCD )
            sigma_relative(0)   = sig_out
            KF_outgoing( 1, 0 ) = kf1
            KF_outgoing( 2, 0 ) = kf2
            if( i_inelastic == 1 )then
                ! Non-zero mass.
                dm1 = m1
                dm2 = m2
                ! Running alpha_s.
                if( alpha_s0 < 1D-10 )then
                    Q2 = s - dm1**2 - dm2**2
                    Nf = 3
                    if( Q2 > amass(4)*amass(4) ) Nf = 4
                    if( Q2 > amass(5)*amass(5) ) Nf = 5
                    alpha_s = func_alpha_s( Nf, Q2, Lambda_QCD2 )
                end if
!       6. Inelastic process: q1 + q1bar -> q2 + q2bar.
                i_proc = 6
                do i_flavor = 1, n_flavor_sample, 1
                    ! KF = 1 : d; KF = 2 : u.
                    i_flavor_out = i_flavor
                    if( i_flavor == 1 ) i_flavor_out = 2
                    if( i_flavor == 2 ) i_flavor_out = 1
                    ! Threshold energy.
                    if( sqrt_s >= 2D0*amass( i_flavor_out ) )then
                        if( ikf1 == i_flavor_out ) cycle
                        dm3    = amass( i_flavor_out )
                        dm4    = amass( i_flavor_out )
                        t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
                        t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
                        t_pi_relative( i_flavor )  = t_min
                        t_0_relative(  i_flavor )  = t_max
                        t_min_relative( i_flavor ) = t_min
                        t_max_relative( i_flavor ) = t_max
                        sig_out = get_sigma( t_min, t_max, t_cut, s,      &
                                             dm1, dm3, alpha_s, K_factor, &
                                             i_proc, i_LOpQCD )
                        sigma_relative( i_flavor ) = sig_out
                        KF_outgoing( 1, i_flavor ) = i_flavor_out * SIGN(1,kf1)
                        KF_outgoing( 2, i_flavor ) = i_flavor_out * SIGN(1,kf2)
                    end if
                end do
!       7. Inelastic process: q1 + q1bar -> g + g.
                i_proc = 7
                i_flavor = 7
                n_flavor_sample = 7
                dm3    = 0D0
                dm4    = 0D0
                t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
                t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
                t_pi_relative( i_flavor )  = t_min
                t_0_relative(  i_flavor )  = t_max
                t_max = t_max - t_cut
                t_min = t_min + t_cut
                t_min_relative( i_flavor ) = t_min
                t_max_relative( i_flavor ) = t_max
                sig_out = get_sigma( t_min, t_max, t_cut, s,      &
                                     dm1, dm2, alpha_s, K_factor, &
                                     i_proc, i_LOpQCD )
                sigma_relative( i_flavor ) = sig_out
                KF_outgoing( 1, i_flavor ) = 21
                KF_outgoing( 2, i_flavor ) = 21
            end if
!           Calculates the total corss section and determines the process to be
!            happened according to the relative cross section.
            sig = 0D0
            do i = 0, n_flavor_sample, 1
                sig = sig + sigma_relative(i)
            end do
            rand_num = PYR(1)
            prob_low = 0D0
            prob_upp = 0D0
            do i = 0, n_flavor_sample, 1
                prob_upp = prob_upp + sigma_relative(i) / sig
                if( rand_num > prob_low .AND. rand_num < prob_upp ) exit
                prob_low = prob_upp
            end do
            i_proc = 5
            if( i == 7 )then
                i_proc = 7
            else if( i > 0 )then
                i_proc = 6
            end if
            kf3 = KF_outgoing( 1, i )
            kf4 = KF_outgoing( 2, i )
            dm3 = amass( kf3 )
            dm4 = amass( kf4 )
            t_min = t_min_relative( i )
            t_max = t_max_relative( i )
            t_pi  = t_pi_relative( i )
            t_0   = t_0_relative( i )
!       8. Elastic process: q1 + q2 -> q1 + q2 / q1    + q2bar -> q1    + q2bar
!                                              / q1bar + q2bar -> q1bar + q2bar
        else if( l_q1q2 .OR. l_q1q2bar .OR. l_q1barq2bar )then
            i_proc = 8
            dm3   = dm1
            dm4   = dm2
            t_min = func_t_pi( s, dm1, dm2, dm3, dm4 )
            t_max = func_t_0(  s, dm1, dm2, dm3, dm4 )
            t_pi  = t_min
            t_0   = t_max
            if( i_LOpQCD > 2 )then
                t_max = t_max - t_cut
            end if
            sig = get_sigma( t_min, t_max, t_cut, s, dm3, dm4, &
                             alpha_s, K_factor, i_proc, i_LOpQCD )
            kf3 = kf1
            kf4 = kf2
        end if
!------------------------   Cross Section Calculating   ------------------------
!-------------------------------------------------------------------------------


        if( sig <= 0D0 ) return
!       0.0389379 (~0.04) is the transformation factor from GeV^-2 to fm^2.
        sig = sig * 0.0389379D0
!       tsmp: the t value sampled.
        tsmp = sample_t( t_min, t_max, t_cut, s, dm1, dm2, dm3, dm4, &
                         i_proc, i_LOpQCD )
!       Symmetric scattering angle.
        ! if( PYR(1) > 0.5D0 ) tsmp = t_min + t_max - tsmp
!       The cosine of the scattering angle theta_s.
        cos_theta_s = 1D0 - 2D0 * ( t_0 - tsmp ) / ( t_0 - t_pi )
!       Equivalently.
!       cos_theta_s = 1D0 + 2D0*s*tsmp / ( xs - am1 - am2 ) / ( xs - am1 + am2 )
        m3 = dm3
        m4 = dm4

        !Lei_debug
        ! write(1100+i_proc,*) s, sig*10D0
        ! if( i_proc /= 1 ) sig = 0D0
        ! adj1(13) = adj1(13) + 1D0
        ! write(*,*) "<< adj1(13)",  adj1(13)   !Lei_debug
        ! call PYFILL( 1, tsmp, 1D0 )
        !Lei_debug

!       Labels heavy quarks realted proceses.
        ikf3 = ABS(kf3)
        ikf4 = ABS(kf4)
        ! Double-heavy, 20+
        if(      ( ikf1 > 3 .AND. ikf2 > 3 .AND. ikf1 /= 21 .AND. ikf2 /= 21 ) &
            .OR. ( ikf3 > 3 .AND. ikf4 > 3 .AND. ikf3 /= 21 .AND. ikf4 /= 21 ) &
            )then
            i_proc = i_proc + 20
        ! Single-heavy, 10+
        else if( ( ikf1 > 3 .AND. ikf1 /= 21 ) &
            .OR. ( ikf2 > 3 .AND. ikf2 /= 21 ) )then
            i_proc = i_proc + 10
        end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function get_sigma( t_min, t_max, t_cut, s, m1, m2, &
                            alpha_s, K_factor, i_proc, i_LOpQCD )
!!      Total cross section of q / qbar /g + q / qbar / g .
!       In this function, we do not distinguish between q and Q.
!       i_proc: internal order number of the process.
!           1: g     + g     -> g     + g
!           2: g     + g     -> q1    + q1bar   (inelastic)
!           3: q     + g     -> q     + g
!              qbar  + g     -> qbar  + g
!           4: q1    + q1    -> q1    + q1
!              q1bar + q1bar -> q1bar + q1bar
!           5: q1    + q1bar -> q1    + q1bar
!           6: q1    + q1bar -> q2    + q2bar   (inelastic)
!           7: q1    + q1bar -> g     + g       (inelastic)
!           8: q1    + q2    -> q1    + q2
!              q1    + q2bar -> q1    + q2bar
!              q1bar + q2bar -> q1bar + q2bar
!#TODO(Lei20241113): Heavy onia gluon dissociation.
!            xx: hOnia + g  -> hQ + hQbar       (dissociation)
!#TODO(Lei20241113): Heavy onia gluon regeneration.
!            xx: hQ + hQbar -> hOnia + g        (regeneration)
        implicit none
        real(kind=8), parameter :: pi=3.141592653589793D0
        real(kind=8), intent(in) :: t_min, t_max, t_cut, s, m1, m2
        real(kind=8), intent(in) :: alpha_s, K_factor
        integer, intent(in) :: i_proc, i_LOpQCD
        real(kind=8) :: sigma
        real(kind=8) :: dME_PF0
        real(kind=8) :: dME_PF1_stu, dME_PF1_sut
        real(kind=8) :: dME_PF1_ust, dME_PF1_uts
        real(kind=8) :: dME_PF2_stu, dME_PF2_ust
        real(kind=8) :: dME_PF3_stu, dME_PF3_tsu
        real(kind=8) :: dME_PF0_REG
        real(kind=8) :: dME_PF1_stu_REG, dME_PF1_sut_REG
        real(kind=8) :: dME_PF1_ust_REG, dME_PF1_uts_REG
        real(kind=8) :: dME_PF2_stu_REG, dME_PF2_ust_REG
        real(kind=8) :: dME_PF3_stu_REG, dME_PF3_tsu_REG
        real(kind=8) :: get_sigma


        get_sigma = 0D0
        sigma = 0D0
        if( t_min > t_max ) return

!       Elastic processes.

!       Small angle and zero-mass approximations for light quarks in
!        elastic parton-parton cross sections with the regulator.
        if( i_LOpQCD == 1 )then
            select case( i_proc )
            case(1)
                sigma = 9D0/2D0 * s**2 * ( 1D0/t_cut - 1D0/(s + t_cut) )
            case(3)
                sigma = 2D0     * s**2 * ( 1D0/t_cut - 1D0/(s + t_cut) )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF3_tsu_REG( s, t_max, MAX( m1, m2 ), t_cut ) &
                          - dME_PF3_tsu_REG( s, t_min, MAX( m1, m2 ), t_cut )
                    sigma = - 8D0 / 3D0 * sigma
                end if
            case( 4 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/(s + t_cut) )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_stu_REG( s, t_max, m1, m1, t_cut ) &
                          - dME_PF1_stu_REG( s, t_min, m1, m1, t_cut ) &
                          + dME_PF1_sut_REG( s, t_max, m1, m1, t_cut ) &
                          - dME_PF1_sut_REG( s, t_min, m1, m1, t_cut ) &
                          + dME_PF2_stu_REG( s, t_max, m1, t_cut ) &
                          - dME_PF2_stu_REG( s, t_min, m1, t_cut )
                    ! Symmetry factor for identical particles in final state.
                    sigma = sigma / 2D0
                end if
            case( 5 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/(s + t_cut) )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_ust_REG( s, t_max, m1, m1, t_cut ) &
                          - dME_PF1_ust_REG( s, t_min, m1, m1, t_cut ) &
                          + dME_PF1_uts_REG( s, t_max, m1, m1, t_cut ) &
                          - dME_PF1_uts_REG( s, t_min, m1, m1, t_cut ) &
                          + dME_PF2_ust_REG( s, t_max, m1, t_cut ) &
                          - dME_PF2_ust_REG( s, t_min, m1, t_cut )
                end if
            case( 8 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/(s + t_cut) )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_stu_REG( s, t_max, m1, m2, t_cut ) &
                          - dME_PF1_stu_REG( s, t_min, m1, m2, t_cut )
                end if
            end select

!       Full LO-pQCD and massive forms for parton-parton cross sections
!        with the regulator.
        else if( i_LOpQCD == 2 )then
            select case( i_proc )
            case(1)
                sigma = dME_PF0_REG( s, t_max, t_cut ) &
                      - dME_PF0_REG( s, t_min, t_cut )
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            case(3)
                sigma = dME_PF3_tsu_REG( s, t_max, MAX( m1, m2 ), t_cut ) &
                      - dME_PF3_tsu_REG( s, t_min, MAX( m1, m2 ), t_cut )
                sigma = - 8D0 / 3D0 * sigma
            case(4)
                sigma = dME_PF1_stu_REG( s, t_max, m1, m1, t_cut ) &
                      - dME_PF1_stu_REG( s, t_min, m1, m1, t_cut ) &
                      + dME_PF1_sut_REG( s, t_max, m1, m1, t_cut ) &
                      - dME_PF1_sut_REG( s, t_min, m1, m1, t_cut ) &
                      + dME_PF2_stu_REG( s, t_max, m1, t_cut ) &
                      - dME_PF2_stu_REG( s, t_min, m1, t_cut )
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            case(5)
                sigma = dME_PF1_ust_REG( s, t_max, m1, m1, t_cut ) &
                      - dME_PF1_ust_REG( s, t_min, m1, m1, t_cut ) &
                      + dME_PF1_uts_REG( s, t_max, m1, m1, t_cut ) &
                      - dME_PF1_uts_REG( s, t_min, m1, m1, t_cut ) &
                      + dME_PF2_ust_REG( s, t_max, m1, t_cut ) &
                      - dME_PF2_ust_REG( s, t_min, m1, t_cut )
            case(8)
                sigma = dME_PF1_stu_REG( s, t_max, m1, m2, t_cut ) &
                      - dME_PF1_stu_REG( s, t_min, m1, m2, t_cut )
            end select

!       Small angle and zero-mass approximations for light quarks in elastic
!        parton-parton cross sections with cutoffs in the integral limits.
        else if( i_LOpQCD == 3 )then
            select case( i_proc )
            case(1)
                sigma = 9D0/2D0 * s**2 * ( 1D0/t_cut - 1D0/s )
            case(3)
                sigma = 2D0     * s**2 * ( 1D0/t_cut - 1D0/s )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF3_tsu( s, t_max, MAX( m1, m2 ) ) &
                          - dME_PF3_tsu( s, t_min, MAX( m1, m2 ) )
                    sigma = - 8D0 / 3D0 * sigma
                end if
            case( 4 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/s )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_stu( s, t_max, m1, m1 ) &
                          - dME_PF1_stu( s, t_min, m1, m1 ) &
                          + dME_PF1_sut( s, t_max, m1, m1 ) &
                          - dME_PF1_sut( s, t_min, m1, m1 ) &
                          + dME_PF2_stu( s, t_max, m1 ) &
                          - dME_PF2_stu( s, t_min, m1 )
                    ! Symmetry factor for identical particles in final state.
                    sigma = sigma / 2D0
                end if
            case( 5 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/s )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_ust( s, t_max, m1, m1 ) &
                          - dME_PF1_ust( s, t_min, m1, m1 ) &
                          + dME_PF1_uts( s, t_max, m1, m1 ) &
                          - dME_PF1_uts( s, t_min, m1, m1 ) &
                          + dME_PF2_ust( s, t_max, m1 ) &
                          - dME_PF2_ust( s, t_min, m1 )
                end if
            case( 8 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/s )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_stu( s, t_max, m1, m2 ) &
                          - dME_PF1_stu( s, t_min, m1, m2 )
                end if
            end select

!       Full LO-pQCD and massive forms for parton-parton cross sections
!        with cutoffs in the integral limits.
        else
            select case( i_proc )
            case(1)
                sigma = dME_PF0( s, t_max ) &
                      - dME_PF0( s, t_min )
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            case(3)
                sigma = dME_PF3_tsu( s, t_max, MAX( m1, m2 ) ) &
                      - dME_PF3_tsu( s, t_min, MAX( m1, m2 ) )
                sigma = - 8D0 / 3D0 * sigma
            case(4)
                sigma = dME_PF1_stu( s, t_max, m1, m1 ) &
                      - dME_PF1_stu( s, t_min, m1, m1 ) &
                      + dME_PF1_sut( s, t_max, m1, m1 ) &
                      - dME_PF1_sut( s, t_min, m1, m1 ) &
                      + dME_PF2_stu( s, t_max, m1 ) &
                      - dME_PF2_stu( s, t_min, m1 )
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            case(5)
                sigma = dME_PF1_ust( s, t_max, m1, m1 ) &
                      - dME_PF1_ust( s, t_min, m1, m1 ) &
                      + dME_PF1_uts( s, t_max, m1, m1 ) &
                      - dME_PF1_uts( s, t_min, m1, m1 ) &
                      + dME_PF2_ust( s, t_max, m1 ) &
                      - dME_PF2_ust( s, t_min, m1 )
            case(8)
                sigma = dME_PF1_stu( s, t_max, m1, m2 ) &
                      - dME_PF1_stu( s, t_min, m1, m2 )
            end select
        end if

!       Inelastic processes.

!       Full LO-pQCD and massive forms for inelastic processes
!        with the regulator.
        if( i_LOpQCD == 1 .OR. i_LOpQCD == 2 )then
            select case( i_proc )
            case(2)
                sigma = dME_PF3_stu_REG( s, t_max, m1, t_cut ) &
                      - dME_PF3_stu_REG( s, t_min, m1, t_cut )
            case(6)
                sigma = dME_PF1_ust_REG( s, t_max, m1, m2, t_cut ) &
                      - dME_PF1_ust_REG( s, t_min, m1, m2, t_cut )
            case(7)
                sigma = dME_PF3_stu_REG( s, t_max, m1, t_cut ) &
                      - dME_PF3_stu_REG( s, t_min, m1, t_cut )
                sigma = 64D0 / 9D0 * sigma
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            end select

!       Full LO-pQCD and massive forms for inelastic processes
!        with cutoffs in the integral limits.
        else
            select case( i_proc )
            case(2)
                sigma = dME_PF3_stu( s, t_max, m1 ) &
                      - dME_PF3_stu( s, t_min, m1 )
            case(6)
                sigma = dME_PF1_ust( s, t_max, m1, m2 ) &
                      - dME_PF1_ust( s, t_min, m1, m2 )
            case(7)
                sigma = dME_PF3_stu( s, t_max, m1 ) &
                      - dME_PF3_stu( s, t_min, m1 )
                sigma = 64D0 / 9D0 * sigma
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            end select
        end if

        if( sigma > 0D0 )  get_sigma = K_factor * sigma * pi * alpha_s**2 &
                                     / (s - (m1 + m2)**2 ) / (s - (m1 - m2)**2 )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function sample_t( t_min, t_max, t_cut, s, m1, m2, m3, m4, &
                           i_proc, i_LOpQCD ) result(t)
!!      Samples t from the differential cross section distribution dsigma/dt.
!       In this function, we do not distinguish between q and Q.
!       i_proc: internal order number of the process.
!           1: g     + g     -> g     + g
!           2: g     + g     -> q1    + q1bar   (inelastic)
!           3: q     + g     -> q     + g
!              qbar  + g     -> qbar  + g
!           4: q1    + q1    -> q1    + q1
!              q1bar + q1bar -> q1bar + q1bar
!           5: q1    + q1bar -> q1    + q1bar
!           6: q1    + q1bar -> q2    + q2bar   (inelastic)
!           7: q1    + q1bar -> g     + g       (inelastic)
!           8: q1    + q2    -> q1    + q2
!              q1    + q2bar -> q1    + q2bar
!              q1bar + q2bar -> q1bar + q2bar
!#TODO(Lei20241113): Heavy onia gluon dissociation.
!            xx: hOnia + g  -> hQ + hQbar       (dissociation)
!#TODO(Lei20241113): Heavy onia gluon regeneration.
!            xx: hQ + hQbar -> hOnia + g        (regeneration)
        implicit none
        integer, parameter :: n_burn_in = 1000
        real(kind=8), intent(in) :: t_min, t_max, t_cut, s, m1, m2, m3, m4
        integer, intent(in) :: i_proc, i_LOpQCD
        ! integer :: n_walk, i
        real(kind=8) :: t, rand_num, rand_num_1, rand_num_2 !, t_propose
        real(kind=8) :: dsigma_dt, dsigma_dt_max
        real(kind=8) :: dsigma_dt_func, PYR
        real(kind=8) :: prob_accept
        ! real(kind=8) :: get_sigma
        ! real(kind=8) :: weight1, weight2, max1, max2, max_seperate, t_separate
        ! real(kind=8) :: norm1, norm2, M_norm


        !Lei_debug
        ! integer :: iHist(100), n_sample, n_fill, i_sample, i_proc, i_LOpQCD
        ! real(kind=8) :: t_min, t_max, t_cut, s, m1, m2, m3, m4
        ! iHist = [ (i, i=1,100,1) ]
        ! i_LOpQCD = 2
        ! i_proc = 1
        ! s = 100D0
        ! t_min = - s
        ! t_max = 0
        ! if( i_LOpQCD == 3 )then
        !     t_max = - t_cut
        ! else if( i_LOpQCD == 4 )then
        !     t_min = - s + t_cut
        !     t_max = - t_cut
        ! end if
        ! ! write(*,*) "t_min, t_max", t_min, t_max
        ! m1 = 0D0
        ! m2 = 0D0
        ! m3 = 0D0
        ! m4 = 0D0
        ! n_sample = 0
        ! n_fill = 1000000
        ! i_sample = 1
        ! i_sample = 2
        !Lei_debug

        t = 0D0

!       Elastic processes wiht the small angle and zero-mass light quarks.
        ! Direct sampling method.
        select case( i_proc )
        case( 1, 3, 4, 5, 8 )
            select case( i_LOpQCD )
            ! Cross sections regularized by the regulator.
            case(1)
                ! With the common form of f(t) = 1/(t - t_cut)^2.
                rand_num = PYR(1)
                t = t_cut * ( 1D0 - ( s + t_cut ) / ( s*rand_num + t_cut ) )
                if( m1+m2 < 1D-1 ) return
            ! Cross sections regularized by cutoffs in the integral limits.
            case(3)
                ! With the common form of f(t) = 1/t^2.
                rand_num = PYR(1)
                t = t_min * t_max / ( t_max + (t_min - t_max)*rand_num )
                if( m1+m2 < 1D-1 ) return
                return
            end select
        end select


        !Lei_debug
        ! return
        ! iQCD2
        ! do while(.true.)
        !     rand_num_1 = PYR(1)
        !     t = t_min + rand_num_1 * ( t_max - t_min )
        !     dsigma_dt     = dsigma_dt_func( s, t, t_cut, m1, m2,    &
        !                                     m3,m4, i_proc, i_LOpQCD )
        !     dsigma_dt_max = MAX(                                    &
        !                 dsigma_dt_func( s, t_min, t_cut, m1, m2,    &
        !                                 m3, m4, i_proc, i_LOpQCD ), &
        !                 dsigma_dt_func( s, t_max, t_cut, m1, m2,    &
        !                                 m3, m4, i_proc, i_LOpQCD ) )
        !     prob_accept = dsigma_dt / dsigma_dt_max
        !     rand_num_2  = PYR(1)
        !     if( rand_num_2 <= prob_accept )then
        !         n_sample = n_sample + 1
        !         call PYFILL(1,t,1D0)
        !         if( n_sample == n_fill )then
        !             call PYDUMP(3,1001,1,iHist(1))
        !             close(1001)
        !             stop
        !         end if
        !     end if
        ! end do
        ! iQCD4
        ! n_walk = 0
        ! do while(.true.)
        !     n_walk = n_walk + 1
        !     rand_num_1  = PYR(1)
        !     t_propose   = t_min + rand_num_1 * ( t_max - t_min )
        !     prob_accept = dsigma_dt_func( s, t_propose, t_cut, &
        !                     m1, m2, m3, m4, i_proc, i_LOpQCD ) &
        !                 / dsigma_dt_func( s, t, t_cut, m1, m2, &
        !                             m3, m4, i_proc, i_LOpQCD )
        !     rand_num_2  = PYR(1)
        !     if( rand_num_2 <= prob_accept )then
        !         t = t_propose
        !         if( n_walk > n_burn_in )then
        !             n_sample = n_sample + 1
        !             call PYFILL(1,t,1D0)
        !             if( n_sample == n_fill )then
        !                 call PYDUMP(3,1001,1,iHist(1))
        !                 close(1001)
        !                 stop
        !             end if
        !         end if
        !     end if
        ! end do
        !Lei_debug


        ! Firstly, tries the rejection sampling "n_burn_in" times.
        ! do i = 1, n_burn_in+1, 1
        ! select case( i_proc )
        ! case( 1 )
        !     t_separate = ( t_min + t_max ) / 2D0
        !     weight1 = get_sigma( t_min, t_separate, t_cut, s, m1, m2, &
        !                          1D0, 1D0, i_proc, i_LOpQCD)
        !     weight2 = get_sigma( t_separate, t_max, t_cut, s, m1, m2, &
        !                          1D0, 1D0, i_proc, i_LOpQCD)
        !     weight1 = weight1 / ( weight1 + weight2 )
        !     weight2 = 1D0 - weight1
        !     max1 =         dsigma_dt_func( s, t_min, t_cut, m1, m2,      &
        !                                    m3, m4, i_proc, i_LOpQCD )
        !     max_seperate = dsigma_dt_func( s, t_separate, t_cut, m1, m2, &
        !                                    m3, m4, i_proc, i_LOpQCD )
        !     max2 =         dsigma_dt_func( s, t_max, t_cut, m1, m2,      &
        !                                    m3, m4, i_proc, i_LOpQCD )
        !     norm1 = 1D0 / ( t_min + s + t_cut )  - 1D0 / ( t_separate +s+t_cut )
        !     norm2 = 1D0 / ( t_separate - t_cut ) - 1D0 / ( t_max - t_cut )
        !     do while(.true.)
        !         rand_num   = PYR(1)
        !         rand_num_1 = PYR(1)
        !         ! Large angle side, i.e. u -> 0. ~ 1/(t + s + t_cut)^2.
        !         if( rand_num <= weight1 )then
        !             t = ( ( t_min + s + t_cut ) &
        !               / ( 1D0 - rand_num_1*(t_min + s + t_cut) ) - s - t_cut ) &
        !               / norm1
        !             dsigma_dt = dsigma_dt_func( s, t, t_cut, m1, m2,    &
        !                                         m3,m4, i_proc, i_LOpQCD )
        !             M_norm = max_seperate &
        !                    * ( t_separate +s +t_cut ) * (t_separate +s +t_cut)
        !             prob_accept = dsigma_dt &
        !                         * ( t + s + t_cut ) * ( t + s + t_cut ) / M_norm
        !         ! Small angle side, i.e. t -> 0. ~ 1/(t - t_cut)^2.
        !         else
        !             t = ( ( t_separate - t_cut ) &
        !               / ( 1D0 - rand_num_1*( t_separate - t_cut ) ) + t_cut ) &
        !               / norm2
        !             dsigma_dt = dsigma_dt_func( s, t, t_cut, m1, m2,    &
        !                                         m3,m4, i_proc, i_LOpQCD )
        !             M_norm = max_seperate &
        !                    * ( t_separate - t_cut ) * ( t_separate - t_cut )
        !             prob_accept = dsigma_dt &
        !                         * ( t - t_cut) * ( t - t_cut ) / M_norm
        !         end if
        !         rand_num_2  = PYR(1)
        !         if( rand_num_2 <= prob_accept )then
        !             n_sample = n_sample + 1
        !             call PYFILL(1,t,1D0)
        !             ! write(*,*) n_sample, t, t_min, t_max, t_separate, max1, max2, max_seperate, M_norm
        !             ! write(*,*) norm1, norm2
        !             if( n_sample == n_fill )then
        !                 call PYDUMP(3,1001,1,iHist(1))
        !                 close(1001)
        !                 stop
        !             end if
        !         end if
        !     end do
        ! end select


!       Full LO-pQCD forms for el. and inel. cross sections.
        ! Firstly, tries the rejection sampling "n_burn_in" times.
        ! do i = 1, n_burn_in+1, 1
        do while(.true.)
            rand_num_1 = PYR(1)
            t = t_min + rand_num_1 * ( t_max - t_min )
            dsigma_dt     = dsigma_dt_func( s, t, t_cut, m1, m2,    &
                                            m3,m4, i_proc, i_LOpQCD )
            dsigma_dt_max = MAX(                                    &
                        dsigma_dt_func( s, t_min, t_cut, m1, m2,    &
                                        m3, m4, i_proc, i_LOpQCD ), &
                        dsigma_dt_func( s, t_max, t_cut, m1, m2,    &
                                        m3, m4, i_proc, i_LOpQCD ) )
            prob_accept = dsigma_dt / dsigma_dt_max
            rand_num_2  = PYR(1)
            if( rand_num_2 <= prob_accept ) exit
        end do
        ! If that fials, uses the Metropolis-Hastings algorithm.
        ! t = t_min + rand_num_1 * ( t_max - t_min )
        ! if( i >= n_burn_in )then
            ! n_walk = 0
            ! do while(.true.)
            !     n_walk = n_walk + 1
            !     rand_num_1  = PYR(1)
            !     t_propose   = t_min + rand_num_1 * ( t_max - t_min )
            !     prob_accept = dsigma_dt_func( s, t_propose, t_cut, &
            !                     m1, m2, m3, m4, i_proc, i_LOpQCD ) &
            !                 / dsigma_dt_func( s, t, t_cut, m1, m2, &
            !                             m3, m4, i_proc, i_LOpQCD )
            !     rand_num_2  = PYR(1)
            !     if( rand_num_2 <= prob_accept )then
            !         t = t_propose
            !         if( n_walk > n_burn_in ) exit
            !     end if
            ! end do
        ! end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dsigma_dt_func( s, t, t_cut, m1, m2, m3, m4, &
                 i_proc, i_LOpQCD )
!!      The differential cross section f(t) = dsigma/dt.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
!       In this function, we do not distinguish between q and Q.
!       i_proc: internal order number of the process.
!           1: g     + g     -> g     + g
!           2: g     + g     -> q1    + q1bar   (inelastic)
!           3: q     + g     -> q     + g
!              qbar  + g     -> qbar  + g
!           4: q1    + q1    -> q1    + q1
!              q1bar + q1bar -> q1bar + q1bar
!           5: q1    + q1bar -> q1    + q1bar
!           6: q1    + q1bar -> q2    + q2bar   (inelastic)
!           7: q1    + q1bar -> g     + g       (inelastic)
!           8: q1    + q2    -> q1    + q2
!              q1    + q2bar -> q1    + q2bar
!              q1bar + q2bar -> q1bar + q2bar
!#TODO(Lei20241113): Heavy onia gluon dissociation.
!            xx: hOnia + g  -> hQ + hQbar       (dissociation)
!#TODO(Lei20241113): Heavy onia gluon regeneration.
!            xx: hQ + hQbar -> hOnia + g        (regeneration)
        implicit none
        integer, intent(in) :: i_proc, i_LOpQCD
        real(kind=8), intent(in) :: s, t, t_cut, m1, m2, m3, m4
        real(kind=8) :: u, dsigma_dt_func
        real(kind=8) :: dME_G0,     dME_G1,     dME_G2,     dME_G3
        real(kind=8) :: dME_G0_REG, dME_G1_REG, dME_G2_REG, dME_G3_REG


        u = m1**2 + m2**2 + m3**2 + m4**4 - s - t
        dsigma_dt_func = 0D0

!       With regulator.
        if( i_LOpQCD == 1 .OR. i_LOpQCD == 2 )then
            select case( i_proc )
            case(1)
                dsigma_dt_func = dME_G0_REG( s, t, u, t_cut )
            case(2)
                dsigma_dt_func = dME_G3_REG( s, t, u, m3, t_cut )
            case(3)
                dsigma_dt_func = - 8D0/3D0 &
                               * dME_G3_REG( t, s, u, MAX(m1,m2), -t_cut )
            case(4)
                dsigma_dt_func = dME_G1_REG( s, t, u, m1,    m1, t_cut ) &
                               + dME_G1_REG( s, u, t, m1,    m1, t_cut ) &
                               + dME_G2_REG( s, t, u, m1, t_cut, t_cut )
            case(5)
                dsigma_dt_func = dME_G1_REG( u, s, t, m1,     m1, -t_cut ) &
                               + dME_G1_REG( u, t, s, m1,     m1,  t_cut ) &
                               + dME_G2_REG( u, s, t, m1, -t_cut,  t_cut )
            case(6)
                dsigma_dt_func = dME_G1_REG( u, s, t, m1, m3, -t_cut )
            case(7)
                dsigma_dt_func = 64D0 / 9D0 * dME_G3_REG( s, t, u, m1, t_cut )
            case(8)
                dsigma_dt_func = dME_G1_REG( s, t, u, m1, m2, t_cut )
            end select
!       Without regulator.
        else
            select case( i_proc )
            case(1)
                dsigma_dt_func = dME_G0( s, t, u )
            case(2)
                dsigma_dt_func = dME_G3( s, t, u, m3 )
            case(3)
                dsigma_dt_func = - 8D0/3D0 * dME_G3( t, s, u, MAX(m1,m2) )
            case(4)
                dsigma_dt_func = dME_G1( s, t, u, m1, m1 ) &
                               + dME_G1( s, u, t, m1, m1 ) &
                               + dME_G2( s, t, u, m1 )
            case(5)
                dsigma_dt_func = dME_G1( u, s, t, m1, m1 ) &
                               + dME_G1( u, t, s, m1, m1 ) &
                               + dME_G2( u, s, t, m1 )
            case(6)
                dsigma_dt_func = dME_G1( u, s, t, m1, m3 )
            case(7)
                dsigma_dt_func = 64D0 / 9D0 * dME_G3( s, t, u, m1 )
            case(8)
                dsigma_dt_func = dME_G1( s, t, u, m1, m2 )
            end select
        end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_alpha_s( Nf, Q2, Lambda_QCD2 )
!!      First-order running strong couping constant alpha_s.
        implicit none
        real(kind=8), parameter :: pi = 3.141592653589793D0
        integer, intent(in) :: Nf
        real(kind=8), intent(in) :: Q2, Lambda_QCD2
        real(kind=8) :: func_alpha_s


        func_alpha_s = 0D0
        if( Q2 <= Lambda_QCD2 ) return
        func_alpha_s = 12D0 * pi / ( 33D0 - 2D0*Nf ) / LOG( Q2 / Lambda_QCD2 )
        ! Upper limit ?
        ! if( func_alpha_s > 1D0 ) func_alpha_s = 1D0


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_mD2( alpha_s, Nf, T, i_statistic )
!!      Debye mass of the gluon within a thermal medium.
!       T: temperature
!       i_statistic: Bose-Einstein and Fermi-Dirac statistics or Boltzmann
!       Cf. Zhe Xu, Phys.Rev.C 71 (2005) 064901.
!           Stephen M.H. Wong, Phys.Rev.C 54 (1996) 2588-2599.
        implicit none
        real(kind=8), parameter :: pi = 3.141592653589793D0
        integer, parameter :: Nc = 3
        real(kind=8), intent(in) :: alpha_s, T
        integer, intent(in) :: Nf, i_statistic
        real(kind=8) :: func_mD2


!       Quantum statistics case, i.e. Bose-Einstein and Fermi-Dirac
!        distributions for gluons and quarks, respectively.
        if( i_statistic == 1 )then
            func_mD2 = 4D0 * pi * alpha_s / 3D0 * ( 1D0*Nc + Nf/2D0 ) * T*T

!       Classical statistics case, i.e. Boltzmann distribution for
!        both gluons and quarks.
        else
            func_mD2 = 8D0 * alpha_s / pi * ( 1D0*Nc + 1D0*Nf ) * T*T
        end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_mq2( alpha_s, T, i_statistic )
!!      Quark medium mass within a thermal medium.
!       T: temperature
!       i_statistic: Bose-Einstein and Fermi-Dirac statistics or Boltzmann
!       Cf. Zhe Xu, Phys.Rev.C 71 (2005) 064901.
!           Stephen M.H. Wong, Phys.Rev.C 54 (1996) 2588-2599.
        implicit none
        real(kind=8), parameter :: pi = 3.141592653589793D0
        integer, parameter :: Nc = 3
        real(kind=8), intent(in) :: alpha_s, T
        integer, intent(in) :: i_statistic
        real(kind=8) :: func_mq2


!       Quantum statistics case, i.e. Bose-Einstein and Fermi-Dirac
!        distributions for gluons and quarks, respectively.
        if( i_statistic == 1 )then
            func_mq2 = pi * alpha_s / 4D0 * ( 1D0*Nc*Nc - 1D0 ) / Nc * T*T

!       Classical statistics case, i.e. Boltzmann distribution for
!        both gluons and quarks.
        else
            func_mq2 = 2D0 * alpha_s / pi * ( 1D0*Nc*Nc - 1D0 ) / Nc * T*T
        end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_t_0( s, m1, m2, m3, m4 )
!!      t_pi < t < t_0: the momentum transfer interval without cutoff.
        implicit none
        real(kind=8), intent(in) :: s, m1, m2, m3, m4
        real(kind=8) :: func_t_0, func_lambda, pow2


        func_t_0 = 1D0 / 4D0 / s * pow2( ( m1*m1 - m2*m2 - m3*m3 + m4*m4 ) )  &
                 - 1D0 / 4D0 / s * pow2( ( func_lambda( s, m1*m1, m2*m2 )     &
                 -                         func_lambda( s, m3*m3, m4*m4 ) ) )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_t_pi( s, m1, m2, m3, m4 )
!!      t_pi < t < t_0: the momentum transfer interval without cutoff.
        implicit none
        real(kind=8), intent(in) :: s, m1, m2, m3, m4
        real(kind=8) :: func_t_pi, func_lambda, pow2


        func_t_pi = 1D0 / 4D0 / s * pow2( ( m1*m1 - m2*m2 - m3*m3 + m4*m4 ) )  &
                  - 1D0 / 4D0 / s * pow2( ( func_lambda( s, m1*m1, m2*m2 )     &
                  +                         func_lambda( s, m3*m3, m4*m4 ) ) )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_lambda( a, b, c )
!!      Kallen (triangle, lambda) function.
        implicit none
        real(kind=8), intent(in) :: a, b, c
        real(kind=8) :: func_lambda


        func_lambda = SQRT( a*a + b*b + c*c - 2D0*a*b - 2D0*a*c - 2D0*b*c )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G0( s, t, u )
!!      The Matrix Element component G0.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
        implicit none
        real(kind=8), intent(in) :: s, t, u
        real(kind=8) :: dME_G0


        dME_G0 = 9D0/2D0 * ( 3D0 - u*t/s**2 - s*u/t**2 - t*s/u**2 )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G0_REG( s, t, u, t_cut )
!!      The Matrix Element component G0 with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, u, t_cut
        real(kind=8) :: dME_G0_REG


        ! dME_G0_REG = 9D0/2D0 * ( 3D0 - u*t/s**2 &
        dME_G0_REG = 9D0/2D0 * ( 3D0 - u*t/(s + t_cut)**2 &
                   - s*u/(t - t_cut)**2 - t*s/(u - t_cut)**2 )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G1( s, t, u, m1, m2 )
!!      The Matrix Element component G1.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
        implicit none
        real(kind=8), intent(in) :: s, t, u, m1, m2
        real(kind=8) :: dME_G1


        dME_G1 = 2D0/9D0 * 2D0/t**2 * ( (s - (m1**2 + m2**2) )**2 &
               + (m1**2 + m2**2 - u)**2 + 2D0*t*(m1**2 + m2**2) )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G1_REG( s, t, u, m1, m2, t_cut )
!!      The Matrix Element component G1(s,t,u).
        implicit none
        real(kind=8), intent(in) :: s, t, u, m1, m2, t_cut
        real(kind=8) :: dME_G1_REG


        dME_G1_REG = 2D0/9D0 * 2D0/(t - t_cut)**2    &
                   * ( (s - (m1**2 + m2**2) )**2     &
                   + (m1**2 + m2**2 - u)**2 + 2D0*t*(m1**2 + m2**2) )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G2( s, t, u, m )
!!      The Matrix Element component G2.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
        implicit none
        real(kind=8), intent(in) :: s, t, u, m
        real(kind=8) :: dME_G2


        dME_G2 = - 2D0/27D0 * 4D0/t/u * (s - 2D0*m**2) * (s - 6D0*m**2)


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G2_REG( s, t, u, m, t_cut, u_cut )
!!      The Matrix Element component G2.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
        implicit none
        real(kind=8), intent(in) :: s, t, u, m, t_cut, u_cut
        real(kind=8) :: dME_G2_REG


        dME_G2_REG = - 2D0/27D0 * 4D0 / (t - t_cut) / (u - u_cut) &
                   * (s - 2D0*m**2) * (s - 6D0*m**2)


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G3( s, t, u, m )
!!      The Matrix Element component G3.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
        implicit none
        real(kind=8), intent(in) :: s, t, u, m
        real(kind=8) :: dME_G3


        dME_G3 = 3D0/16D0 * 4D0 * (m**2 - t) * (m**2 - u) / s**2          &
               + 1D0/12D0 * ( 2D0 * (m**2 - u) * (m**2 - t)               &
                - 4D0 * m**2 * (m**2 + t) ) / (t - m**2)**2               &
               + 1D0/12D0 * ( 2D0 * (m**2 - t) * (m**2 - u)               &
                - 4D0 * m**2 * (m**2 + u) ) / (u - m**2)**2               &
               - 1D0/96D0 * 4D0 * m**2 * (s - 4D0*m**2)                   &
                / (m**2 - t) / (m**2 - u)                                 &
               + 3D0/32D0 * 4D0 * ( t * (s + t) - m**4 ) / s / (m**2 - t) &
               + 3D0/32D0 * 4D0 * ( u * (s + u) - m**4 ) / s / (m**2 - u)


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G3_REG( s, t, u, m, s_cut )
!!      The Matrix Element component G3.
        implicit none
        real(kind=8), intent(in) :: s, t, u, m, s_cut
        real(kind=8) :: dME_G3_REG


        dME_G3_REG = 3D0/16D0 * 4D0 * (m**2 - t) * (m**2 - u) / (s + s_cut)**2 &
                   + 1D0/12D0 * ( 2D0 * (m**2 - u) * (m**2 - t)                &
                    - 4D0 * m**2 * (m**2 + t) ) / (t - m**2)**2                &
                   + 1D0/12D0 * ( 2D0 * (m**2 - t) * (m**2 - u)                &
                    - 4D0 * m**2 * (m**2 + u) ) / (u - m**2)**2                &
                   - 1D0/96D0 * 4D0 * m**2 * (s - 4D0*m**2)                    &
                     / (m**2 - t) / (m**2 - u)                                 &
                   + 3D0/32D0 * 4D0 * ( t * (s + t) - m**4 ) &
                    / ( s + s_cut ) / ( m**2 - t ) &
                   + 3D0/32D0 * 4D0 * ( u * (s + u) - m**4 ) &
                    / ( s + s_cut ) / ( m**2 - u )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF0( s, t )
!!      The primitive function of the Matrix Element component G0.
        implicit none
        real(kind=8), intent(in) :: s, t
        real(kind=8) :: dME_PF0


        dME_PF0 = 3D0*t + t**2/2D0/s + t**3/3D0/s**2 - s**2/t &
                + s*LOG( ABS( t / (s + t) ) ) - s**2/(s + t)
        dME_PF0 = 9D0 / 2D0 * dME_PF0


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF0_REG( s, t, t_cut )
!!      The primitive function of the Matrix Element component G0 w/ regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, t_cut
        real(kind=8) :: dME_PF0_REG


        dME_PF0_REG = 3D0*t + s*t**2/2D0/(s + t_cut)**2             &
                    + t**3/3D0/(s + t_cut)**2                       &
                    - s**2/(t - t_cut)                              &
                    + s*LOG( ABS( (t - t_cut) / (s + t + t_cut) ) ) &
                    - s*t_cut/(t - t_cut) - (s**2 + s*t_cut)/(s + t + t_cut)
        dME_PF0_REG = 9D0 / 2D0 * dME_PF0_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_stu( s, t, m1, m2 )
!!      The primitive function of the Matrix Element component G1(s,t,u).
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2
        real(kind=8) :: dME_PF1_stu, M


        M = m1**2 + m2**2
        dME_PF1_stu = - 2D0*(s - M)**2 / t &
                    + 2D0 * S * LOG( ABS(t) ) + t
        dME_PF1_stu = 4D0 / 9D0 * dME_PF1_stu


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_stu_REG( s, t, m1, m2, t_cut )
!!      The primitive function of the Matrix Element component
!!       G1(s,t,u) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2, t_cut
        real(kind=8) :: dME_PF1_stu_REG, M


        M = m1**2 + m2**2
        dME_PF1_stu_REG = - 2D0*(s - M)**2 / (t - t_cut)  &
                        - 2D0*(s - M)*t_cut / (t - t_cut) &
                        - t_cut**2 / (t - t_cut)          &
                        + (t - t_cut)                     &
                        + (2D0*s + 2D0*t_cut) * LOG( ABS(t - t_cut) )
        dME_PF1_stu_REG = 4D0 / 9D0 * dME_PF1_stu_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_sut( s, t, m1, m2 )
!!      The primitive function of the Matrix Element component G1(s,u,t).
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2
        real(kind=8) :: dME_PF1_sut, M


        M = m1**2 + m2**2
        dME_PF1_sut = 2D0*(s - M)**2 / (2D0*M - s - t) &
                    - (2D0*M - s - t)                  &
                    - 2D0*s*LOG( ABS(2D0*M - s - t) )
        dME_PF1_sut = 4D0 / 9D0 * dME_PF1_sut


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_sut_REG( s, t, m1, m2, t_cut )
!!      The primitive function of the Matrix Element component
!!       G1(s,u,t) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2, t_cut
        real(kind=8) :: dME_PF1_sut_REG, M


        M = m1**2 + m2**2
        dME_PF1_sut_REG = ( (s - M)**2 + (s - M + t_cut)**2 ) &
                          / (2D0*M - s - t - t_cut)           &
                        - (2D0*M - s - t - t_cut)             &
                        - 2D0*(s + t_cut)*LOG( ABS(2D0*M - s - t - t_cut) )
        dME_PF1_sut_REG = 4D0 / 9D0 * dME_PF1_sut_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_ust( s, t, m1, m2 )
!!      The primitive function of the Matrix Element component G1(u,s,t).
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2
        real(kind=8) :: dME_PF1_ust, M


        M = m1**2 + m2**2
        dME_PF1_ust = 2D0 * (t - M)**3 / 3D0 / s**2 &
                    + t**2 / s + t
        dME_PF1_ust = 4D0 / 9D0 * dME_PF1_ust


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_ust_REG( s, t, m1, m2, s_cut )
!!      The primitive function of the Matrix Element component
!!       G1(s,u,t) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2, s_cut
        real(kind=8) :: dME_PF1_ust_REG, M


        M = m1**2 + m2**2
        dME_PF1_ust_REG = 2D0 * (t - M)**3 / 3D0 / (s + s_cut)**2 &
                        + t**2 / (s + s_cut) + s**2 * t / (s + s_cut)**2
        dME_PF1_ust_REG = 4D0 / 9D0 * dME_PF1_ust_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_uts( s, t, m1, m2 )
!!      The primitive function of the Matrix Element component G1(u,t,s).
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2
        real(kind=8) :: dME_PF1_uts, M


        M = m1**2 + m2**2
        dME_PF1_uts = -2D0*(M - s)**2 / t + 2D0*s*LOG( ABS(t) ) + t
        dME_PF1_uts = 4D0 / 9D0 * dME_PF1_uts


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_uts_REG( s, t, m1, m2, t_cut )
!!      The primitive function of the Matrix Element component
!!       G1(u,t,s) with the regulator
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2, t_cut
        real(kind=8) :: dME_PF1_uts_REG, M


        M = m1**2 + m2**2
        dME_PF1_uts_REG = -( (M - s)**2 + (M - s - t_cut)**2 ) / (t - t_cut) &
                        + 2D0*(s + t_cut)*LOG( ABS(t - t_cut) )              &
                        + (t - t_cut)
        dME_PF1_uts_REG = 4D0 / 9D0 * dME_PF1_uts_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF2_stu( s, t, m1 )
!!      The primitive function of the Matrix Element component G2(s,t,u).
        implicit none
        real(kind=8), intent(in) :: s, t, m1
        real(kind=8) :: dME_PF2_stu, M


        M = m1**2 + m1**2
        dME_PF2_stu = (s - M) * (s - 3D0*M) &
                    * LOG( ABS( t / (2D0*M - s - t) ) ) / (2D0*M - s)
        dME_PF2_stu = - 8D0 / 27D0 * dME_PF2_stu


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF2_stu_REG( s, t, m1, t_cut )
!!      The primitive function of the Matrix Element component
!!       G2(s,t,u) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, t_cut
        real(kind=8) :: dME_PF2_stu_REG, M


        M = m1**2 + m1**2
        dME_PF2_stu_REG = (s - M) * (s - 3D0*M) &
                        * LOG( ABS( (t - t_cut) / (2D0*M - s - t - t_cut) ) ) &
                        / (2D0*M - s - 2D0*t_cut)
        dME_PF2_stu_REG = - 8D0 / 27D0 * dME_PF2_stu_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF2_ust( s, t, m1 )
!!      The primitive function of the Matrix Element component G2(u,s,t).
        implicit none
        real(kind=8), intent(in) :: s, t, m1
        real(kind=8) :: dME_PF2_ust, M


        M = m1**2 + m1**2
        dME_PF2_ust = 2D0*t + t**2/2D0/s &
                    + (s - M**2/s) * LOG( ABS(t) )
        dME_PF2_ust = - 8D0 / 27D0 * dME_PF2_ust


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF2_ust_REG( s, t, m1, t_cut )
!!      The primitive function of the Matrix Element component
!!       G2(u,s,t) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, t_cut
        real(kind=8) :: dME_PF2_ust_REG, M


        M = m1**2 + m1**2
        dME_PF2_ust_REG = ( 2D0*(s + t_cut)*(t - t_cut)          &
                        + 1D0/2D0*(t - t_cut)**2                 &
                        + (s**2 - M**2 + 2D0*s*t_cut + t_cut**2) &
                        * LOG( ABS(t - t_cut) ) )                &
                        / (s + t_cut)
        dME_PF2_ust_REG = - 8D0 / 27D0 * dME_PF2_ust_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF3_stu( s, t, m1 )
!!      The primitive function of the Matrix Element component G3(s,t,u).
        implicit none
        real(kind=8), intent(in) :: s, t, m1
        real(kind=8) :: dME_PF3_stu, T1, T2, T3, T4, T5, T6, M


        M  = m1**2
        T1 = 4D0 * ( M*t**2 - s*t**2/2D0 &
           - t**3/3D0 - M**2*t + M*s*t ) / s**2
        T2 = 2D0 * ( - t - s*LOG( ABS(M - t) ) ) &
           - 4D0 * M * ( 2D0*M/(M - t) + LOG( ABS(M - t) ) )
        T3 = 2D0 * ( s*LOG( ABS(s + t - M) ) - t ) &
           - 4D0 * M * ( 2D0*M / (M - s - t) - LOG( ABS(M - s - t) ) )
        T4 = 4D0 * M*(s - 4D0*M)/s * LOG( ABS( (s + t - M) / (M - t) ) )
        T5 = 4D0 * ( s*M - s*t + 2D0*M**2 - 2D0*M*t &
           - (M - t)**2 / 2D0 - s*M*LOG( ABS(M - t) ) ) / s
        T6 = 4D0 * ( (s + t - M)**2 / 2D0 - s**2 - s*t - s*M - 2D0*M*t &
           + 2D0*M**2 + s*M*LOG( ABS(s + t - M) ) ) / s
        ! T1 = 0D0
        dME_PF3_stu = 3D0/16D0*T1 + 1D0/12D0*T2 &
                    + 1D0/12D0*T3 - 1D0/96D0*T4 &
                    + 3D0/32D0*T5 + 3D0/32D0*T6


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF3_stu_REG( s, t, m1, s_cut )
!!      The primitive function of the Matrix Element component
!!       G3(s,t,u) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, s_cut
        real(kind=8) :: dME_PF3_stu_REG, T1, T2, T3, T4, T5, T6, M


        M  = m1**2
        T1 = 4D0 * ( M*t**2 - s*t**2/2D0 &
           - t**3/3D0 - M**2*t + M*s*t ) / (s + s_cut)**2
        T2 = 2D0 * ( - t - s*LOG( ABS(M - t) ) ) &
           - 4D0 * M * ( 2D0*M/(M - t) + LOG( ABS(M - t) ) )
        T3 = 2D0 * ( s*LOG( ABS(s + t - M) ) - t ) &
           - 4D0 * M * ( 2D0*M / (M - s - t) - LOG( ABS(M - s - t) ) )
        T4 = 4D0 * M*(s - 4D0*M)/s * LOG( ABS( (s + t - M) / (M - t) ) )
        T5 = 4D0 * ( s*M - s*t + 2D0*M**2 - 2D0*M*t &
           - (M - t)**2 / 2D0 - s*M*LOG( ABS(M - t) ) ) / (s + s_cut)
        T6 = 4D0 * ( (s + t - M)**2 / 2D0 - s**2 - s*t - s*M - 2D0*M*t &
           + 2D0*M**2 + s*M*LOG( ABS(s + t - M) ) ) / (s + s_cut)
        dME_PF3_stu_REG = 3D0/16D0*T1 + 1D0/12D0*T2 &
                        + 1D0/12D0*T3 - 1D0/96D0*T4 &
                        + 3D0/32D0*T5 + 3D0/32D0*T6


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF3_tsu( s, t, m1 )
!!      The primitive function of the Matrix Element component G3(t,s,u).
        implicit none
        real(kind=8), intent(in) :: s, t, m1
        real(kind=8) :: dME_PF3_tsu, T1, T2, T3, T4, T5, T6, M


        M  = m1**2
        T1 = 4D0 * (M - s) * ( M/t - s/t + LOG( ABS(t) ) )
        T2 = 2D0 * ( s*t + t**2/2D0 - M*t ) / (M - s) &
           - 4D0 * M * (M + s) / (s - M)**2 * t
        T3 = 2D0 * (M - s) * LOG( ABS(s + t - M) ) &
           - 4D0 * M * ( 2D0*M / (M - s - t) - LOG( ABS(M - s - t) ) )
        T4 = 4D0 * M / (M - s) * ( s + t - M &
           - s*LOG( ABS(s + t - M) ) - 3D0*M*LOG( ABS(s + t - M) ) )
        T5 = 4D0 * ( s*t + (s**2 - M**2) * LOG( ABS(t) ) ) / (M - s)
        T6 = 4D0 * ( (3D0*M**2 - 4D0*M*s + s**2) / (s - M) &
           * LOG( ABS( t / ( s + t - M ) ) ) &
           - (2D0*M - s) * LOG( ABS(s + t - M) ) )
        dME_PF3_tsu = 3D0/16D0*T1 + 1D0/12D0*T2 &
                    + 1D0/12D0*T3 - 1D0/96D0*T4 &
                    + 3D0/32D0*T5 + 3D0/32D0*T6


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF3_tsu_REG( s, t, m1, t_cut )
!!      The primitive function of the Matrix Element component G3(t,s,u).
        implicit none
        real(kind=8), intent(in) :: s, t, m1, t_cut
        real(kind=8) :: dME_PF3_tsu_REG, T1, T2, T3, T4, T5, T6, A, B, C, D, M


        M  = m1**2
        T1 = 4D0 * (M - s) * ( s/(t_cut - t) + t_cut/(t_cut - t) &
           + LOG( ABS(t - t_cut) ) + M/(t - t_cut) )
        T2 = 2D0 * ( s*t + t**2/2D0 - M*t ) / (M - s) &
           - 4D0 * M * (M + s) / (s - M)**2 * t
        T3 = 2D0 * (M - s) * LOG( ABS(s + t - M) ) &
           - 4D0 * M * ( 2D0*M / (M - s - t) - LOG( ABS(M - s - t) ) )
        T4 = 4D0 * M / (M - s) * ( s + t - M &
           - s*LOG( ABS(s + t - M) ) - 3D0*M*LOG( ABS(s + t - M) ) )
        T5 = 4D0 * ( -M**2*LOG( ABS(t - t_cut) ) &
           + s * ( (s + t_cut)*LOG( ABS(t - t_cut) ) + (t - t_cut) ) ) &
           / (M - s)
         A = ( 3D0*M**2 - 4D0*M*s + s**2 ) / ( s + t_cut - M )
         B = - A
         C = - t_cut * ( 2D0*M - s ) / ( s + t_cut - M )
         D = - ( 2D0*M - s ) - C
        T6 = 4D0 * ( (A + C) * LOG( ABS(t - t_cut) ) &
           + (B + D) * LOG( ABS(s + t - M) ) )
        dME_PF3_tsu_REG = 3D0/16D0*T1 + 1D0/12D0*T2 &
                        + 1D0/12D0*T3 - 1D0/96D0*T4 &
                        + 3D0/32D0*T5 + 3D0/32D0*T6


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine collis( icp )
!!      Performs parton-parton collision & updates particle list
!       icp: the current icp-th collision pair
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (MSCA=20000)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/papr_p/core,xs,xu,xt,sm,as,dta,xa,sl0,tl0,qa, &
         ea,sqq,sqg,sgg,pa(3),pip(10,msca),mtime,kfk,nsca,kpip(msca)
        common/shower_1/ KSH(6,5), PSH(6,5), VSH(6,5)
        common/syspar_p/rsig1,pio,tcut
        common/trs/ntrs,nontrs,ktrs(kszj,5),ptrs(kszj,5),vtrs(kszj,5)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/sa25/i_inel_proc,i_time_shower,i_diquark,ipad25,para1_1,para1_2
        common/sa26/ndiq(kszj),npt(kszj),ifcom(kszj),idi,idio
        common/sa28/nstr,nstra(kszj),nstrv(kszj),nstr0, &
         nstr1,nstr1a(kszj),nstr1v(kszj)
        common/work7/reac(100),crose(100)
        common/ctllist_p/nreac(100),nrel
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        common/scatt/ ic, jc, n0
        common/e_cons1/ pss(kszj,5), ff(kszj)
!       lc(1,i), lc(2,i): line numbers in the particle list of i-th coll. pair
!       lc(3,i), lc(4,i): flavor codes of scattered particles of i-th coll. pair
!       lc(5,i): the internal process number
!       lc(6,i): event type. =1, elastic scattering; =2, inelastic; =3, shower.
!       tc: the collision time of the colliding pair
!       sig_tsmp_Cthetas_p(1,i): the total cross section of the colliding pair
!       sig_tsmp_Cthetas_p(2,i): the sampled mom. transfer "tsmp" of the pair
!       sig_tsmp_Cthetas_p(3,i): the scattering angle COS(theta_s) of the pair
!       sig_tsmp_Cthetas_p(4,i): the momentum modulus of th pair in the CM frame
!       ic, jc: line number of colliding particles
!       ik1, ik2: line number of the colliding pair in parton list
!       kf3, kf4: kf codes of the collided pair after the collision
!       ifcom(i): line number (in 'sbe') of first component of i-th diquark
!       nreac(i): statistics of the # of successful i-th collision
!       nrel: statistics of the # of collision blocked
!       nsca: number of particles after collision
!       pip(1-3,*): momentum of particle after collision
!       pip(4,*): energy of particle after collision
!       pip(5,*): virtuality of particle after collision
!       pip(6,*): x value of particle after collision
!       kpip(*): flavor code of particle after collision
        integer, intent(in) :: icp
        dimension pi(4), pj(4), pi_o(4), pj_o(4), ps(4)
        dimension px(4), py(4), pij(4), b(3)


!       Gets contents of the current collision pair.
        ic     = lc(1,icp)
        jc     = lc(2,icp)
        kf3    = lc(3,icp)
        kf4    = lc(4,icp)
        i_proc = lc(5,icp)
        i_type = lc(6,icp)
        time   = tc(1,icp)
        sig         = sig_tsmp_Cthetas_p(1,icp)
        tsmp        = sig_tsmp_Cthetas_p(2,icp)
        cos_theta_s = sig_tsmp_Cthetas_p(3,icp)
        p_modulus   = sig_tsmp_Cthetas_p(4,icp)
        ik1 = ic
        ik2 = jc

!       The hadronization model.
        i_had_model = INT( adj1(12) )
!       kf1 and kf2: kf code of the colliding pair before interaction
        kf1 = K(ik1,2)
        kf2 = K(ik2,2)
!       Selection of the form of differential cross section in the LO-pQCD.
        i_LOpQCD = INT( adj1(20) )


!-------------------------------------------------------------------------------
!----------------------------   Paticle Collision   ----------------------------
!       Number of scattered (+ showered) particles.
        nsca  = 0
        pi(4) = P(ik1,4)
        pj(4) = P(ik2,4)
        do i=1,3,1
            pi(i) = P(ik1,i)
            pj(i) = P(ik2,i)
            b(i)  = ( pi(i) + pj(i) ) / ( pi(4) + pj(4) )
        end do
        pi_o = pi
        pj_o = pj
        pij  = pi + pj
!       Invariant mass squared of the colliding pair.
        eiej2 = pij(4)**2 - pij(1)**2 - pij(2)**2 - pij(3)**2
!       Boosts to the current CMS of the colliding pair.
        call lorntz( 0, b, pi, pj )
        dm1 = P(ik1,5)
        dm2 = P(ik2,5)
        dm3 = amass( kf3 )
        dm4 = amass( kf4 )
        am1 = dm1
        am2 = dm2
        am3 = dm3
        am4 = dm4
!       Small angle and zero-mass of light quarks approximations for el. scat.
        if( i_LOpQCD == 1 .OR. i_LOpQCD == 3 )then
            i_proc1 = MOD( i_proc, 10 )
            select case( i_proc1 )
            case( 1, 3, 4, 5, 8 )
                if( ABS(kf1) < 4 ) am1 = 0D0
                if( ABS(kf2) < 4 ) am2 = 0D0
                if( ABS(kf3) < 4 ) am3 = 0D0
                if( ABS(kf4) < 4 ) am4 = 0D0
            end select
        end if
!       Mandelstam variables.
        xs  = eiej2
        xt  = tsmp
        xu  = am1**2 + am2**2 + am3**2 + am4**2 - xs - xt
!       Squareed momentum transfer (Q^2) to be used in the time-like shower.
!       pT^2, same as that in PYTHIA ( Q^2_{hard} ).
        qc2 = xt * xu / xs
!       qc2 = 4D0 * qc2
!       qc2 = MIN( -xt, -xu )

!       Calculates the direction cosines of momentum of one colliding
!        particle after scattering relative to the momentum
!        of the corresponding particle before scattering.
        Ctheta_s = cos_theta_s
!       Random 2*pi azimuthal angle.
        phi_s = 2D0 * pio * PYR(1)

        if( ABS(Ctheta_s) > 1D0 )then
            Ctheta_s = Ctheta_s / ABS( Ctheta_s ) - 1D-10
        end if
        Stheta_s = SQRT( 1D0 - Ctheta_s**2 )
        Cphi_s   = COS( phi_s )
        Sphi_s   = SIN( phi_s )

!       pi and pj, as input into "rotate", are four momentum of the colliding
!        pair (before scattering and rotation), as output are four
!        momentum of that pair after scattering and rotation
!       p_modulus: momentum modulus of pi or pj, both are equal in their CMS.
        call rotate( Ctheta_s, Stheta_s, Cphi_s, Sphi_s, p_modulus, pi, pj )
        pi(4) = SQRT( pi(1)*pi(1) + pi(2)*pi(2) + pi(3)*pi(3) + dm3*dm3 )
        pj(4) = SQRT( pj(1)*pj(1) + pj(2)*pj(2) + pj(3)*pj(3) + dm4*dm4 )

!       Boosts back.
        call lorntz( 1, b, pi, pj )

!       Final state of two parton scattering (nsca=2).
!       Also the initial state of the time-like branching.
        nsca = 2
        kpip(1) = kf3
        kpip(2) = kf4
        do i=1,4,1
            pip(i,1) = pi(i)
            pip(i,2) = pj(i)
        end do
        pip(5,1) = qc2
        pip(5,2) = qc2
        pip(6,1) = 1D0
        pip(6,2) = 1D0

!       Updates particle list (note: line numbers of the colliding
!        pair after interaction are the same as ones before interaction).
        ! ik1
        K(ik1,2) = kpip(1)
        do k1=1,4,1
            P(ik1,k1) = pip(k1,1)
        end do
!       For the inelastic processes.
        if( kf3 /= kf1 ) P(ik1,5) = dm3
        ! ik2
        K(ik2,2) = kpip(2)
        do k1=1,4,1
            P(ik2,k1) = pip(k1,2)
        end do
!       For the inelastic processes.
        if( kf4 /= kf2 ) P(ik2,5) = dm4

!       Counts cross sections.
        reac(i_proc)  = reac(i_proc)  + 1D0
        crose(i_proc) = crose(i_proc) + sig
        nreac(i_proc) = nreac(i_proc) + 1
!----------------------------   Paticle Collision   ----------------------------
!-------------------------------------------------------------------------------


!       Without the time-like shower.
        if( i_time_shower == 0 ) return

        ! Backup for shower.
        ! 1
        KSH(1,2) = K(ik1,2)
        KSH(1,3) = ik1
        do i=1,5,1
            PSH(1,i) = P(ik1,i)
        end do
        ! 2
        KSH(2,2) = K(ik2,2)
        KSH(2,3) = ik2
        do i=1,5,1
            PSH(2,i) = P(ik2,i)
        end do


!-------------------------------------------------------------------------------
!----------------------------   Time-like Shower   -----------------------------
!       n_call_tili: statistics of the times calling "ti_li1".
        n_call_tili = 0
!       tl0: the virtuality cutoff of the time-like branching, i. e. Mu0**2.
        tl0 = adj1(24)

!       Prepars the time-like branching.
        nsca0 = 0
        kmsca = nsca
        do i = nsca0+1, kmsca, 1
            kfk = kpip(i)
            if( ABS(kfk) > 21 ) cycle
            if( pip(5,i) > 4D0*tl0 )then
!       Phase space for branching vanishes at a virtuality of 4*tl0.
!        cf. Eq.25 in B.R. Webber, Ann. Rev. Nucl. Part. Sci. 36 (86) 253.
                ea = pip(4,i)
                qa = pip(5,i)
                do j=1,3
                    pa(j) = pip(j,i)
                end do
                xa  = pip(6,i)
!       Time-like branching for the final state partons & the induced
!        partons formed in the time like branchings.
                    call ti_li1( n_call_tili, i )
                    n_call_tili = n_call_tili + 1
!-------update of i---------------------
                    ! kpip(i) = kfk
                    ! do j=1,3
                    !     pip(j,i) = pa(j)
                    ! enddo
                    ! pip(4,i) = ea
                    ! pip(5,i) = qa
                    ! pip(6,i) = xa
!-------finished------------------------
            end if
        end do
        nsca0 = kmsca
!----------------------------   Time-like Shower   -----------------------------
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!----------------------------   Time-like Shower   -----------------------------
!       Proceeds for the time-like shower.

!       Throws away the old collision pairs.
        call update_ctlm( icp )

!       n_call_tili: statistics of the times calling "ti_li1".
        n_call_tili = 0
!       tl0: the virtuality cutoff of the time-like branching, i. e. Mu0**2.
        tl0 = adj1(24)

!0700   Performing time-like branching ---------------------------------
        nsca0 = 0
! 300   continue
        kmsca = nsca
        do i = nsca0+1, kmsca, 1
            if( pip(5,i) > 4D0*tl0 )then
!       Phase space for branching vanishes at a virtuality of 4*tl0.
!        cf. Eq.25 in B.R. Webber, Ann. Rev. Nucl. Part. Sci. 36 (86) 253.
                ea = pip(4,i)
                qa = pip(5,i)
                do j=1,3
                    pa(j) = pip(j,i)
                end do
                kfk = kpip(i)
                xa  = pip(6,i)
!       Time-like branching for the final state partons & the induced
!        partons formed in the time like branchings.
                    call ti_li1( n_call_tili, i )
                    n_call_tili = n_call_tili + 1
!-------update of i---------------------
                    ! kpip(i) = kfk
                    ! do j=1,3
                    !     pip(j,i) = pa(j)
                    ! enddo
                    ! pip(4,i) = ea
                    ! pip(5,i) = qa
                    ! pip(6,i) = xa
!-------finished------------------------
            end if
        end do
        nsca0 = kmsca
!       if( msca0 /= nsca ) goto 300
!0700   Performing time-like branching finished-------------------------

!0800   "nsca" to "PYJETS" ---------------------------------------------
        if( nsca == 2 )then
!       No partons were showered off. Updates the collision time list.
            call update_ctln
            return
        else if( nsca > 2 )then
!       Branching: g -> gg, g -> qqbar, & q-> qg (qbar -> qbarg)
!       If we consider only g -> gg & no g -> qqbar, it is possible to
!        construct chains of kf3 -> kf3 + kf5, kf4 -> kf4 + kf6
!       if nsca=3, makes "nsca" in order of kf3 (or kf4), kf5, kf6
!       if nsca=4, makes "nsca" in order of kf3, kf4, kf5, kf6
            if( nsca == 3 )then
                do i1=1,nsca,1
                    if( kpip(i1) == kf3 .or. kpip(i1) == kf4 )then
!       Moves particle i1 to the first position "nsca".
                        kpip1   = kpip(1)
                        pip1    = pip(1,1)
                        pip2    = pip(2,1)
                        pip3    = pip(3,1)
                        pip4    = pip(4,1)
                        pip5    = pip(5,1)
                        pip6    = pip(6,1)
                        kpip(1) = kpip(i1)
                        do i2=1,6,1
                            pip(i2,1) = pip(i2,i1)
                        end do
                        kpip(i1)  = kpip1
                        pip(1,i1) = pip1
                        pip(2,i1) = pip2
                        pip(3,i1) = pip3
                        pip(4,i1) = pip4
                        pip(5,i1) = pip5
                        pip(6,i1) = pip6
                    end if
                end do
            else if( nsca == 4 )then
                do i1=1,nsca,1
                    if( kpip(i1) == kf3 )then
!       Moves particle i1 to the first position "nsca".
                        kpip1   = kpip(1)
                        pip1    = pip(1,1)
                        pip2    = pip(2,1)
                        pip3    = pip(3,1)
                        pip4    = pip(4,1)
                        pip5    = pip(5,1)
                        pip6    = pip(6,1)
                        kpip(1) = kpip(i1)
                        do i2=1,6,1
                            pip(i2,1) = pip(i2,i1)
                        end do
                        kpip(i1)  = kpip1
                        pip(1,i1) = pip1
                        pip(2,i1) = pip2
                        pip(3,i1) = pip3
                        pip(4,i1) = pip4
                        pip(5,i1) = pip5
                        pip(6,i1) = pip6
                    end if
                end do
                do i3=2,nsca,1
                    if( kpip(i3) == kf4 )then
!       Moves particle i3 to the second position in "nsca".
                        kpip1   = kpip(2)
                        pip1    = pip(1,2)
                        pip2    = pip(2,2)
                        pip3    = pip(3,2)
                        pip4    = pip(4,2)
                        pip5    = pip(5,2)
                        pip6    = pip(6,2)
                        kpip(2) = kpip(i3)
                        do i2=1,6,1
                            pip(i2,2) = pip(i2,i3)
                        end do
                        kpip(i3)  = kpip1
                        pip(1,i3) = pip1
                        pip(2,i3) = pip2
                        pip(3,i3) = pip3
                        pip(4,i3) = pip4
                        pip(5,i3) = pip5
                        pip(6,i3) = pip6
                    end if
                end do
            end if
            do i=1,3,1
                ! Three position of ik1.
                px(i) = V(ik1,i)
                ! Three position of ik2.
                py(i) = V(ik2,i)
            end do
!       Induced parton is distributed randomly in between colliding pair.
            n0 = N
            do i=3,nsca,1
                rl = PYR(1)
                n0 = n0 + 1
                do k1=1,3
                    V(n0,k1) = px(k1) * rl + py(k1) * ( 1D0 - rl )
                end do
                V(n0,4) = time
                V(n0,5) = 0D0
                ! It will be give below.
                ! taup(n0) = time
                ishp(n0) = 1
            end do

!--------- keeps four momentum conservation ----------------------
            do i2=1,4
                ! pi_o (pj_o) four momentum of colliding pair (ik1 & ik2)
                ps(i2) = pi_o(i2) + pj_o(i2)
            end do
            do i1=1,nsca,1
                do i2=1,4,1
                    pss(i1,i2) = pip(i2,i1)
                enddo
                pss(i1,5) = P(i1,5)
                if( nsca > 2 ) pss(i1,5) = amass( kpip(i1) )
            end do
!       Adjusts four momentum conservation by iteration,no more than
!        5000 iterations.
            call cconse( ps, 1, nsca, 1 )
!       Total momentum of partons after scattering (including induced partons
!        in time-like branching) is normalized to the total momentum of
!        scattering pair.
            do i1=1,nsca,1
                do i2=1,4,1
                    pip(i2,i1) = pss(i1,i2)
                end do
            end do
!       Gives the formation time to the showerd partons.
            n0 = N
            do i1=3,nsca,1
                n0 = n0 + 1
                E   = pip(4,i1)
                pT2 = pip(1,i1)**2 + pip(2,i1)**2
                taup(n0) = time
                !#TODO(Lei20241202)
                ! taup(n0) = time + E / pT2
            end do
!---------- keep four momentum conservation finished -------------

!       A part (first & second) of "nsca" to "PYJETS", i.e. feedbacks
!        first & second partons in "nsca" to ik1 & ik2.
            ! ik1
            K(ik1,2) = kpip(1)
            if( pip(4,1) < 1D-10 ) pip(4,1) = 1D-10
            do k1=1,3,1
                P(ik1,k1) = pip(k1,1)
            end do
            P(ik1,4) = pip(4,1)
            P(ik1,5) = amass( kpip(1) )
            ! ik2
            K(ik2,2) = kpip(2)
            if( pip(4,2) < 1D-10 ) pip(4,2) = 1D-10
            do k1=1,3,1
                P(ik2,k1) = pip(k1,2)
            end do
            P(ik2,4) = pip(4,2)
            P(ik2,5) = amass( kpip(2) )
!0800   finished--------------------------------------------------------

!       Deals with third & fourth partons in "nsca".
            if(nsca == 3)then
                if( i_had_model == 0 .OR. i_had_model == 3 )then
!       Moves third induced partons to "trs" for sfm.
                    ntrs = ntrs + 1
                    ktrs(ntrs,1) = 2
                    ktrs(ntrs,2) = kpip(3)
                    ktrs(ntrs,3) = 0
                    ktrs(ntrs,4) = 0
                    ktrs(ntrs,5) = 0
                    do j=1,4,1
                        ptrs(ntrs,j) = pip(j,3)
                        vtrs(ntrs,j) = V(n0,j)
                    end do
                    ptrs(ntrs,5) = 0D0
                    vtrs(ntrs,5) = 0D0
                elseif( i_had_model /= 0 .AND. i_had_model /= 3 )then
!       Adds this parton (gluon) into collision list for coal.
                    N = N + 1
                    K(N,1) = 2
                    K(N,2) = kpip(3)
                    K(N,3) = 0
                    K(N,4) = 0
                    K(N,5) = 0
                    do j=1,4,1
                        P(N,j) = pip(j,3)
                    end do
                    P(N,5) = 0D0
                end if
            else if(nsca == 4)then
!       Composess third & fourth partons in "nsca" as a string & moves it
!        to the end of "PYJETS".
                ! 3
                N = N + 1
                if( pip(4,3) < 1D-10 ) pip(4,3) = 1D-10
                do k1=1,3
                    P(N,k1) = pip(k1,3)
                end do
                P(N,4) = pip(4,3)
                P(N,5) = amass( kpip(3) )
                K(N,1) = 2
                K(N,2) = kpip(3)
                K(N,3) = 0
                K(N,4) = 0
                K(N,5) = 0
                ! 4
                N = N + 1
                if( pip(4,4) < 1D-10 ) pip(4,4) = 1D-10
                do k1=1,3,1
                    P(N,k1) = pip(k1,4)
                end do
                P(N,4) = pip(4,4)
                P(N,5) = amass( kpip(4) )
                K(N,1) = 1
                K(N,2) = kpip(4)
                K(N,3) = 0
                K(N,4) = 0
                K(N,5) = 0
            end if
        end if

!       Copies the new string composed of the final state parton pair
!        (ik1 & ik2) to the end of "PYJETS".
        if( ik2 > ik1 )then
            K(ik1,1) = 2
            K(ik2,1) = 1
!       Copies ik1 to N+1.
            call coend(ik1)
            icnew = N
!       Copies ik2 to N+1.
            call coend(ik2)
            jcnew = N
        else
            K(ik1,1) = 1
            K(ik2,1) = 2
!       Copies ik2 to N+1.
            call coend(ik2)
            jcnew = N
!       Copies ik1 to N+1.
            call coend(ik1)
            icnew = N
        end if
        nstr1 = nstr1 + 1
        nstr1a(nstr1) = N - 1
        nstr1v(nstr1) = N

!       Handles the scattering parton pair (scattered parton pair has been
!        treated as a new string above).

!       2: g + g -> q + qbar
        if( i_proc == 2 )then
!       Throws away scattering parton pair (ik1 & ik2).
            if(ik1 > ik2)then
                call parmov( ik1, ik1, ik2, i_proc )
                call parmov( ik2, ik1, ik2, i_proc )
            else if(ik1 < ik2)then
                call parmov( ik2, ik1, ik2, i_proc )
                call parmov( ik1, ik1, ik2, i_proc )
            endif
            ! New partons are above n-2 w/ or w/o branching.
            n00 = N - 2
            call updpli_p( n00 )
            return
        end if

!       6: q1 + q1bar -> q2 + q2bar, 7: q + qbar -> g + g
        if( i_proc == 6 .OR. i_proc == 7 )then
!       hadronization with coalescence
            if( i_had_model /= 0 .AND. i_had_model /= 3 )then
!       throws away scattering parton pair (ik1 & ik2)
                if(ik1 > ik2)then
                    call parmov( ik1, ik1, ik2, i_proc )
                    call parmov( ik2, ik1, ik2, i_proc )
                endif
                if(ik1 < ik2)then
                    call parmov( ik2, ik1, ik2, i_proc )
                    call parmov( ik1, ik1, ik2, i_proc )
                endif
                ! New partons are above n-2 w/ or w/o branching.
                n00 = N - 2
                call updpli_p( n00 )
                return
!       hadronization with string fragmentation.
            elseif( i_had_model == 0 .OR. i_had_model == 3 )then
!       is ik1 (ik2) a component of string ?
                call adjst( ik1, ik1str, ik1sa, ik1sv )
!       ik1str: oredr number of string to which ik1 belongs, equal 0 otherwise
!       ik1sa: line number of first component of above string, equal 0 otherwise
!       ik1sv: line number of last component of above string, equal 0 otherwise
                call adjst( ik2, ik2str, ik2sa, ik2sv )
!       is ik1 a component of diquark ?
                call adjdi( ik1, idi1, iway1 )
                ifcom1 = 0
                npt1   = 0
                if(iway1 == 2) npt1 = npt(idi1)
!       npt1: line # (in 'sbe') of partner of idi1-th diquark
                if(iway1 == 1) ifcom1 = ifcom(idi1)
!       ifcom1: line # (in 'sbe') of first component of idi1-th diquark
!       does ik2 is a component of diquark ?
                call adjdi( ik2, idi2, iway2 )
                ifcom2 = 0
                npt2 = 0
                if(iway2 == 2) npt2   = npt(idi2)
                if(iway2 == 1) ifcom2 = ifcom(idi2)
!       ifcom1: line # of first component of idi1-th diquark
!       Removes diquarks (to which ik1 & ik2 belong) from diquark list.
!       It doesn't affect other lists.
                if( (idi1 /= 0 .and. idi2 /= 0) .and. idi1 == idi2 ) &
                    call diqmov(idi1)
                if( (idi1 /= 0 .and. idi2 /= 0) .and. idi1 > idi2 )then
                    call diqmov(idi1)
                    call diqmov(idi2)
                end if
                if( (idi1 /= 0 .and. idi2 /= 0) .and. idi1 < idi2 )then
                    call diqmov(idi2)
                    call diqmov(idi1)
                end if
                if(idi1 /= 0 .and. idi2 == 0) call diqmov(idi1)
                if(idi1 == 0 .and. idi2 /= 0) call diqmov(idi2)
!       Moves components of string out & updates lists.
                if(ik1str == 0 .and. ik2str == 0)then
!       Moves ik1 & ik2 out as well as updates lists.
                    if(ik1 > ik2)then
                        call parmov( ik1, ik1, ik2, i_proc )
                        call parmov( ik2, ik1, ik2, i_proc )
                    end if
                    if(ik1 < ik2)then
                        call parmov( ik2, ik1, ik2, i_proc )
                        call parmov( ik1, ik1, ik2, i_proc )
                    end if
                    n00 = N - 2
                    call updpli_p( n00 )
                    return
                end if
                if(ik1str == 0 .and. ik2str /= 0)then
                    if(ik1 < ik2)then
                        call strmov( ik2str, ik2sa, ik2sv, ik1, ik2, i_proc )
                        call parmov( ik1, ik1, ik2, i_proc )
                    end if
                    if(ik1 > ik2)then
                        call parmov( ik1, ik1, ik2, i_proc )
                        call strmov( ik2str, ik2sa, ik2sv, ik1, ik2, i_proc )
                    end if
                    n00 = N - ( ik2sv - ik2sa + 1 ) - 1
                    call updpli_p( n00 )
                    return
                end if
                if(ik1str /= 0 .and. ik2str == 0)then
                    if(ik1 < ik2)then
                        call parmov( ik2, ik1, ik2, i_proc )
                        call strmov( ik1str, ik1sa, ik1sv, ik1, ik2, i_proc )
                    else
                        call strmov( ik1str, ik1sa, ik1sv, ik1, ik2, i_proc )
                        call parmov( ik2, ik1, ik2, i_proc )
                    endif
                    n00 = N - ( ik1sv - ik1sa + 1 ) - 1
                    call updpli_p( n00 )
                    return
                end if
!       Proceeds for both of ik1str & ik2str are not equal to zero.
                if(ik1str == ik2str)then
                    call strmov( ik1str, ik1sa, ik1sv, ik1, ik2, i_proc )
                    n00 = N - ( ik1sv - ik1sa + 1 )
                    call updpli_p( n00 )
                end if
                if(ik1str > ik2str)then
                    call strmov( ik1str, ik1sa, ik1sv, ik1, ik2, i_proc )
                    call strmov( ik2str, ik2sa, ik2sv, ik1, ik2, i_proc )
                    n00 = N - ( ik1sv - ik1sa + 1 ) - ( ik2sv - ik2sa + 1 )
                    call updpli_p( n00 )
                end if
                if(ik1str < ik2str)then
                    call strmov( ik2str, ik2sa, ik2sv, ik1, ik2, i_proc )
                    call strmov( ik1str, ik1sa, ik1sv, ik1, ik2, i_proc )
                    n00 = N - ( ik1sv - ik1sa + 1 ) - ( ik2sv - ik2sa + 1 )
                    call updpli_p( n00 )
                end if
                return
            end if
        end if
!----------------------------   Time-like Shower   -----------------------------
!-------------------------------------------------------------------------------


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine ti_li1( n_call_tili, ij )
!!      Performs the time-like branching (along main chain until tl0=Mu0**2
!!       for the branching parton).
!       n_call_tili: statistics the times calling 'ti_li1'
!       ij: order # of spliting parton
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (MSCA=20000)
        PARAMETER (pio=3.141592653589793D0)
        common/papr_p/core,xs,xu,xt,sm,as,dta,xa,sl0,tl0,qa, &
         ea,sqq,sqg,sgg,pa(3),pip(10,msca),mtime,kfk,nsca,kpip(msca)
        common/shower_1/ KSH(6,5), PSH(6,5), VSH(6,5)
        common/sa24/adj1(40),nnstop,non24,zstop
        dimension pa1(4),pa2(4),p00(4)
        integer :: n_sample_q


        i_stop = 0
!       i_stop=0 means the branching will go on, i_stop=1 means "stop".
        do i=1,4,1
            pa1(i) = 0D0
            pa2(i) = 0D0
            p00(i) = 0D0
        end do
        n_sample_q = 0
10      continue
        n_sample_q = n_sample_q + 1

!       Varieties in time-like branchs are all the mass squared (m**2)
!        and m**2 = E**2 - |p|**2 > 0
!------------------------------------------------------------------
!       Branching:   A  -> A1   +   A2
!       Vituality:  qqb    qqa     qqc
!       ID No.   :  kfk    kf1     kf2
!        x       :  xa x'  xa1 x   xa2
!       Momentum :  pa     pa1     pa2
!
!       Equation : qqb = qqa/z + qqc/(1-z) + pT2 / ( z*(1 - z) ) (*)
!-------------------------------------------------------------------
!       Above Eq. is deduced from the relation of mass and energy (Einstein
!        relation) and four-momentum conservation, cf. note made in
!        Webber's paper.
!-------------------------------------------------------------------
!       A2 will be always the radiated one.

!-------Flavors of A1, A2 will be decided in the ensuing module.
        q2max = qa
        kf0   = kfk
        call suda_ti( i_fail, kf0, i_split, q2max, q2, z )
!       Time-like Sudakov factor, input: kf0 & q2max, output: others.
!       i_fail: =1, failed in sampling q2 and z
!               =0, succeeded in sampling q2 and z
        if( i_fail == 1 ) then
            ! Forcibly stops.
            qa = tl0
            return
        end if
        qqb = q2
        if( n_call_tili > 2 ) qqb = qa
        zab = z
!       "zab" is the sampled ratio of x/x' in this branching.
        if( ABS(zab) <= 1D-10 ) return
        xa1 = xa*zab
        xa2 = xa - xa1
        kf1 = 0
        kf2 = 0
!       Initial state is quark.
        if( kfk /= 21 ) then
!       Branching q -> qg.
            kf1 = kfk
            kf2 = 21
!       Initial state is gluon.
        else
!       Branching g -> gg.
            if( i_split == 2 ) then
                kf1 = 21
                kf2 = 21
            else
!       Considers g->gg only.
                goto 200
!       Branching g -> qqbar
!       ea = pip(4,ij): energy of spliting particle ij
!       pip(1-3,ij): three momentum of spliting particle ij
                eaa = ea*ea - pip(1,ij)*pip(1,ij) - pip(2,ij)*pip(2,ij) &
                    - pip(3,ij)*pip(3,ij)
                if( eaa <= 0D0 ) goto 200
                ! Invariant mass of the initial state gluon.
                eaa = SQRT(eaa)
                if( eaa < ( 2D0*amass(2) ) ) goto 200
                call break_f(eaa,kff,amq)
                kf1 =  kff
                kf2 = -kf1
200         end if
        end if
!----------Finished.

        q2max = zab*( qqb - tl0 / (1D0 - zab) )
!       The max. value of A1's virtuality will be decided by the eqution(*)
!         with qqc=tl0 & pT2=0 .
        call suda_ti( i_fail, kf1, i_split, q2max, q2, z )
        if( i_fail == 1 ) then
!       qqa=tl0 means the forward evolution has finished, thus set i_stop=1.
            qqa = tl0
            i_stop = 1
        else
            qqa = q2
        end if
        q2max = ( qqb - qqa/zab ) * ( 1D0 - zab )
!       The max. value of A2's virtuality will be decided by the Eq(*) w/ pT2=0.
        call suda_ti( i_fail, kf2, i_split, q2max, q2, z )
        if( i_fail == 1 ) then
            qqc = tl0
        else
            qqc = q2
        endif
        pT2 = ( qqb - qqa/zab - qqc/(1D0 - zab) ) * ( 1D0 - zab ) * zab
!       pT2: the square of transverse momentum decided by the Eq.(*)

!0130   The following block gives the momentum of A1, A2 i.e. pa1, pa2
!        according to the definition of z & pT.
        pa02 = pa(1)**2 + pa(2)**2 + pa(3)**2
        ! Rotates to the frame in which the momentum of A is along +Z direction.
        ! Third momentum of A. pT = 0.
        paz  = SQRT(pa02)
        ! Third momentum of A1.
        pz1  = zab*paz
        ! Third momentum of A2.
        pz2  = ( 1D0 - zab ) * paz
        ! Momentum modulus of A1.
        ppa1 = SQRT( pz1**2 + pT2 )
        ! Momentum modulus of A2.
        ppa2 = SQRT( pz2**2 + pT2 )
        ppt  = SQRT(pT2)

!       Imposes the dead-cone effect.
        theta_split = ATAN( ppt / pz2 )
        E_mother    = PSH(ij,4)
        dm_mother   = PSH(ij,5)
        theta_dead  = dm_mother / E_mother
!       Suppresses the mass-related small angle emmision. Resamples.
        if( n_sample_q > 1000 ) return
        if( theta_split < theta_dead ) goto 10

!       Calculates the formation time of a radiated gluon (quark?)
!        in the Lab frame for the LPM effect.
        E_split1 = SQRT( pz1**2 + pT2 + amass(kf1) )
        t_form1  = 2D0*E_split1 / ( pT2 + zab*dm_mother**2 )
        E_split2 = SQRT( pz2**2 + pT2 + amass(kf2) )
        t_form2  = 2D0*E_split2 / ( pT2 + ( 1D0 - zab )*dm_mother**2 )

        sctas1 = ppt / ppa1
        cctas1 = pz1 / ppa1
        phi = 2D0*pio*PYR(1)
        cfis = COS(phi)
        sfis = SIN(phi)
!       Direction cosines of A1 relative to momentum of A
!       Note : transverse momentum of A1 & A2 must be in the same plane
!        in order to keep the transverse momentum conservation.
        sctas2 = -ppt/ppa2
        cctas2 =  pz2/ppa2
        do i=1,3,1
            pa1(i) = pa(i)
            pa2(i) = pa(i)
        end do
        call rotate( cctas1, sctas1, cfis, sfis, ppa1, pa1, p00 )
!       Originally, momentum of A1 is relative to the direction of
!        the momentum of A
!       After rotation, momentum of A1 is relative to the CMS where A
!        belongs to.
        call rotate( cctas2, sctas2, cfis, sfis, ppa2, pa2, p00 )
!0130   Finished---------------------------------------

!0140   Newly induced parton (A2) should be added to particle list.
        nsca = nsca + 1
        kpip(nsca) = kf2
        pip(4,nsca) = amass(kf2)**2
        do i=1,3
            pip(4,nsca) = pip(4,nsca) + pa2(i)**2
            pip(i,nsca) = pa2(i)
        end do
        pip(4,nsca) = SQRT( pip(4,nsca) )
!       pip(5,nsca) = qa2
        pip(5,nsca) = qqc
        pip(6,nsca) = xa2
!0140   Finished--------------------------------------------

!-------Updates five very important values of KFK, QA, XA, PA, EA (time-like)
!        & then perform the forward evolution along main chain further.
        kfk = kf1
        qa  = qqa
        do i=1,3
            pa(i) = pa1(i)
        enddo
        ea = SQRT( pa(1)**2 + pa(2)**2 + pa(3)**2 + amass(kfk)**2 )
        xa = xa1

!       Records infomation of showered partons.
        ! Mother or the first showerd parton.
        KSH(ij,2) = kf1
        do i=1,3,1
            PSH(ij,i) = pa1(i)
        end do
        PSH(ij,4) = ea
        PSH(ij,5) = amass(kf1)
        VSH(ij,5) = t_form1
        ! The second showered parton.
        KSH( ij+3 ,2 ) = kf2
        KSH( ij+3 ,3 ) = KSH(ij,3)
        do i=1,3,1
            PSH( ij+3, i ) = pa2(i)
        end do
        PSH( ij+3 , 4 ) = pip(4,nsca)
        PSH( ij+3 , 5 ) = amass(kf2)
        VSH( ij+3 , 5 ) = t_form2

!       if( n_sample_q > 10 ) return
        ! Once branching only.
        if( n_sample_q >= 1 ) return
        if( i_stop == 0 ) goto 10
!       i_stop=1 means the evolution "stop" while i_stop=0 "continue"
!       The time-like branching of A2 parton will consider later, i.e.
!        in "collis" after statement numbered 300.


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine suda_ti( i_fail, kf0, i_split, q2max, q2, z )
!!      Time-like Sudakov factor.
!!      Performs time-like forward branching, one step only.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (MSCA=20000)
        PARAMETER (pio=3.141592653589793D0)
        common/papr_p/core,xs,xu,xt,sm,as,dta,xa,sl0,tl0,qa, &
         ea,sqq,sqg,sgg,pa(3),pip(10,msca),mtime,kfk,nsca,kpip(msca)
        common/sa24/adj1(40),nnstop,non24,zstop
        real(kind=8) :: Lamda_QCD2
!----------------------------------------------------------------------
!       process    : A  ->   A1  +  A2
!       virtuality : qqb     qqa    qqc
!
!       equation: qqb = qqa / z + qqc / (1 - z) + pT2 / (1 - z) / z
!----------------------------------------------------------------------


!       yint1, 2: the integral functions of part of the splitting functions
        ! q -> q + g, cf. Sa' note p.44.
        yint1(x) = -LOG( 1D0 - x )
        ! g -> g + g
        yint2(x) =  LOG( 1D0 / (1D0 - x) - 1D0 )
        Lamda_QCD2 = adj1(25)**2
        i_fail = 1
        i_split = 0
        z = 0D0
        ! Cf. W' paper p.264 Eq.25.
        if( q2max <= 4*tl0 ) return
        ikf = ABS(kf0)
        tmax = LOG( q2max / Lamda_QCD2 )
        ! Cf. W' paper p.264
        tmin = LOG( 4D0*tl0 / Lamda_QCD2 )
!       tmax, tmin: bounds of t = dlog( Q**2/lamda ) in the Sudakov form factor
        ! cf. W' paper p.263
        zmin = tl0/q2max
        ! cf. W' paper p.263
        zmax = 1D0 - tl0/q2max
!       Approximated values of the bounds of z.
        if( zmax <= zmin ) return
        ! q -> qg
        if( ikf /= 21 )then
            ! cf. Sa' note p. 43
            ciqg = 4D0/3D0 * ( ( 2D0*yint1(zmax) - zmax - 0.5D0*zmax**2 ) &
                             - ( 2D0*yint1(zmin) - zmin - 0.5D0*zmin**2 ) )
!       ciqg: the integration of splitting function for q -> qg
            cisum = ciqg
            i_split = 1
        ! g ->
        else
            cigg = ( 6D0*yint2(zmax) - 12D0*zmax + 3D0*zmax**2 - 2D0*zmax**3 ) &
                 - ( 6D0*yint2(zmin) - 12D0*zmin + 3D0*zmin**2 - 2D0*zmin**3 )
!       cigg: the integration of the splitting function for g->gg
            ciqq = 1D0/2D0 * ( ( zmax - zmax**2 + 2D0/3D0*zmax**3 ) &
                             - ( zmin - zmin**2 + 2D0/3D0*zmin**3 ) )
!       ciqq: the integration of the splitting function for g->qqbar
            cisum = cigg + ciqq
            aaa = PYR(1)
            ! g -> gg
            if( aaa <= ( cigg / cisum ) ) then
                i_split = 2
                cisum = cigg
            ! g -> qqbar
            else
                i_split = 3
                cisum = ciqq
            end if
        end if
        ce = 9D0/2D0/cisum
        n_sample_z = 0
100     continue
        aaa = PYR(1)
        ! Running alpha_s.
!       tt = tmax * aaa**ce
        ! Constant alpha_s.
        tt = tmax + 2D0*pio / as / cisum * LOG(aaa)
!       tt: the t value sampled from time-like Sudakov factor
        if( tt <= tmin ) return
        q2 = Lamda_QCD2*EXP(tt)
!       Lamda_QCD2: Lambda_s**2
        rzmax = 0.5D0 + 0.5D0 * SQRT( 1D0 - 4D0*tl0/q2 )
!       cf. Eq.24 in B. R. Webber, Ann. Rev. Nucl. Part. Sci. 36(1986)253
!        that journal is simplified as W' elsewhere
        rzmin = 1D0 - rzmax
!       rzmax, rzmin: exact value of zmax & zmin corresponding to the t sampled.
!0170   Samples z. -------------------------------------------------------------
200     continue
        aaa = PYR(1)
        bbb = PYR(1)
        ! q -> qg
        if( i_split == 1 )then
            zsp = aaa*yint1(zmax) + (1D0 - aaa)*yint1(zmin)
            zsp = 1D0 - EXP(-zsp)
            ! The acceptance probability.
            ratio = (1D0 + zsp*zsp) / 2D0
            ! if(ratio <= bbb) goto 200
            ! Same as the previous statement.
            if( bbb > ratio ) goto 200
        ! g -> gg
        else if( i_split == 2 )then
            zsp = aaa*yint2(zmax) + (1D0 - aaa)*yint2(zmin)
            zsp = 1D0 - 1D0 / (1D0 + EXP(zsp) )
            ratio = ( 1D0 - zsp*(1D0 - zsp) )**2
            ! if(ratio <= bbb) goto 200
            ! Same as the previous statement.
            if( bbb > ratio ) goto 200
        ! g -> qqbar
        else
            zsp = aaa*zmax + (1D0 - aaa)*zmin
            ratio = 1D0 - 2D0*zsp + 2D0*zsp**2
            ! if(ratio <= bbb) goto 200
            ! Same as the previous statement.
            if(bbb > ratio) goto 200
        end if
!0170   Sampling z finished. ---------------------------------------------------
        if( zsp > rzmax .OR. zsp < rzmin ) then
!       If the z sampled falls out of the exact region, reject the sampled t.
!       Let tmax equal to the t sampled, and go back to resample.
            n_sample_z = n_sample_z + 1
!       Forcibly stops the time-like branching.
            if( n_sample_z > 1000 )then
                i_fail = 1
                return
            end if
            tmax = tt
            goto 100
        end if
        z = zsp
        i_fail = 0


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine cconse(ps,npl,np,nstep)
!!      Adjusts four momentum conservation by iteration,no more than
!!       5000 iterations
!       pp : four momentum of particles
!       ps : above four momenta should be conserved to ps
!       npl : line # of the first particle
!       np : line # of last particle
!       nstep : interval of the step
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000)
        common/e_cons1/ pp(kszj,5), ff(kszj)
        dimension ps(4),pxyz(3),arp(3)


        dep  = 1D-5
        ps4  = ps(4)
        pxyz = 0D0
        jj = 0
100     es = 0D0
        do i=npl,np,nstep
            es = es + pp(i,4)
        end do
        fr = es / ps4
!       f(dabs(1.-fr)  <=  dep) return
        do i=npl,np,nstep
            amas  = pp(i,5)
            amas2 = amas*amas
            ppm   = pp(i,4)
            ppf   = ppm / fr
            ff(i) = SQRT( ABS( ppf*ppf - amas2 ) / ( ppm*ppm - amas2 ) )
            do j=1,3,1
                ppp = ff(i) * pp(i,j)
                pp(i,j) = ppp
                pxyz(j) = pxyz(j) + ppp
            end do
        end do
        do i=1,3,1
            arp(i) = ABS( 1D0 - pxyz(i) / ps(i) )
        end do
        if( ABS(1D0 - fr) <= dep .AND. arp(1) <= dep .AND. arp(2) <= dep &
            .AND. arp(3) <= dep ) return
        do i=1,3,1
            pxyz(i) = pxyz(i) - ps(i)
            pxyz(i) = pxyz(i) / ( DBLE(np-npl) / DBLE(nstep) + 1D0 )
        end do
        do i=npl,np,nstep
            do j=1,3,1
                pp(i,j) = pp(i,j) - pxyz(j)
            end do
            pp5  = pp(i,5)
            pp52 = pp5 * pp5
            pp(i,4) = SQRT( pp52 + pp(i,1)**2 + pp(i,2)**2 + pp(i,3)**2 )
        end do
        jj = jj + 1
        if(jj < 5000) goto 100


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function fap_qq(zz)
!!      Splitting function for g-->q*q(-)
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)


        fap_qq = ( zz**2 + (1D0 - zz)**2 ) / 2D0


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function fap_gg(zz)
!!      Splitting function for g-->g*g.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)


        fap_gg = 6D0 * ( 1D0 - zz * ( 1D0 - zz ) )**2 / ( zz * ( 1D0 - zz ) )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function fap_qg(zz)
!!      Splitting function for q-->q*g.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)


        fap_qg = 4D0 * ( 1D0 + zz**2 ) / 3D0 / ( 1D0 - zz )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function amass(kf)
!!      Mass of parton.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        common/sa18/i_deex,n_deex_step,i_pT_coal,i_pT_endpoint,a_FF,aPS_c,aPS_b
        common/sa24/adj1(40),nnstop,non24,zstop


!       Hadronization model.
        i_had_model = INT( adj1(12) )
        i_mass = 0
!       For the debug, uses the current algebra mass.
        if( i_had_model == 1 .AND. i_deex > 99 ) i_mass = 3
        amass = 0D0
        r     = 0D0
        k1    = ABS(kf)
        select case( k1 )
        case( 21 )
        case( 2 )
            r = 0.33D0
            if( i_mass == 3 ) r = 0.0056D0
        case( 1 )
            r = 0.33D0
            if( i_mass == 3 ) r = 0.0099D0
        case( 3 )
            r = 0.5D0
            if( i_mass == 3 ) r = 0.199D0
        case( 4 )
            r = 1.5D0
            if( i_mass == 3 ) r = 1.23D0
        case( 5 )
            r = 4.8D0
            if( i_mass == 3 ) r = 4.17D0
        case( 6 )
            r = 175D0
            if( i_mass == 3 ) r = 165D0
        end select
        amass = r


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine coend(ii)
!!      Copies the parton ii to the end of "PYJETS".
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa26/ndiq(kszj),npt(kszj),ifcom(kszj),idi,idio


        N = N + 1
        do i3=1,5,1
            K(N,i3) = K(ii,i3)
            P(N,i3) = P(ii,i3)
            V(N,i3) = V(ii,i3)
        end do
        taup(N) = taup(ii)
        ishp(N) = ishp(ii)
!       ndiq(N) = ndiq(ii)
!       do j1=1,idi,1
!           kk = ifcom(j1)
!           jj = npt(j1)
!           if( kk == ii ) ifcom(j1) = N
!           if( jj == ii ) npt(j1)   = N
!       end do


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine adjst(ik1,ik1str,ik1sa,ik1sv)
!!      Finds order number etc. of ik1 string.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000)
        common/sbe/nbe,nonbe,kbe(kszj,5),pbe(kszj,5),vbe(kszj,5)
        common/sa28/nstr,nstra(kszj),nstrv(kszj),nstr0, &
         nstr1,nstr1a(kszj),nstr1v(kszj)


        ik1str = 0   ! if ik1 not belong to string
        ik1sa  = 0   ! if ik1 not belong to string
        ik1sv  = 0   ! if ik1 not belong to string
        if(ik1 <= nbe)then   ! nbe is given before call break
            do i1 = 1,nstr0,1   ! nstr0 is given after call break
                i1a = nstr1a(i1)   ! i1a: line number of 'A' of i1-th string
                i1v = nstr1v(i1)   ! i1v: line number of 'V' of i1-th string
                if(ik1 >= i1a .and. ik1 < i1v+1)then
                    ik1str = i1    ! order number of string to which ik1 belongs
                    ik1sa  = i1a   ! line number of 'A' of above string
                    ik1sv  = i1v   ! line number of 'V' of above string
                    return
                end if
            end do
        end if
        if(ik1 > nbe .and. nstr1 > nstr0)then
            do i1 = nstr0+1, nstr1, 1
                i1a = nstr1a(i1)
                i1v = nstr1v(i1)
                if( ik1 >= i1a .and. ik1 < i1v+1 )then
                    ik1str = i1
                    ik1sa  = i1a
                    ik1sv  = i1v
                    return
                end if
            end do
        end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine adjdi(ii,idi1,idway)
!!      Is ii component of a diquark.
!       idway=0: if ii is not first component or it's partner of a diquark
!       idway=1: if ii is first component of a diquark
!       idway=2: if ii is partner of a diquark
!       idi1: order # of diquark to which ii belong
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sbe/nbe,nonbe,kbe(kszj,5),pbe(kszj,5),vbe(kszj,5)
        common/sa26/ndiq(kszj),npt(kszj),ifcom(kszj),idi,idio
!       ifcom(i): line number (in 'sbe') of first component of i-th diquark


        idway = 0
        idi1  = 0
        do i1=1,idi,1
            i2 = ifcom(i1)
            if(ii == i2)then
!       ii is first component of idi1-th diquark
                idi1 = i1
                idway = 1
                return
            end if
        end do

        do i1=1,idi,1
            ii1=npt(i1)   ! line number of partner of i1-th diquark
            if(ii == ii1)then   ! ii is partner of i1-th diquark
!       ii is partner of idi1-th diquark
                idi1=i1   ! ii is partner of idi1-th diquark
                idway=2
                return
            end if
        end do


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine diqmov(idi1)
!!      Removes diquark idi1 out of list.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000)
        common/sbe/nbe,nonbe,kbe(kszj,5),pbe(kszj,5),vbe(kszj,5)
        common/sa26/ndiq(kszj),npt(kszj),ifcom(kszj),idi,idio


        if(idi > 0)then
            do i1=1,nbe
                if(ndiq(i1) == idi1) ndiq(i1) = 0
!       ndiq(j): = 0 if j is quark (antiquark)
!                = order # of diquark if j is diquark (anti-diquark)
!             j: line number in 'sbe'
            end do
            if(idi1 == idi) idi = idi - 1
            if(idi1 < idi)then
!       updates diquark list
                do i1 = idi1+1, idi, 1
                    i2 = i1 - 1
                    ifcom(i2) = ifcom(i1)
!       ifcom(idi): line number of first component of idi-th diquark
                    npt(i2)=npt(i1)
                end do
            idi = idi - 1
            end if
        end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine parmov(ii,ik1,ik2, i_proc )
!!      If "ii" is not equal to ik1 (ik2) and i_proc=6 or 7,
!!        moves parton "ii" to "trs" as well as updates lists.
!!      Otherwise, throws away parton ii (ii is equal to ik1 (ik2) ),
!!       i.e. updates lists
!       Updates "PYJETS" one step downward from ii+1 to N.
!       Updates "sbe" one step downward from ii+1 to nbe if ii <= nbe.
!       Updates diquark list/
!       Updates string list.
!       Updates collision time list.
!       ik1 & ik2 are the line number (in parton list) of the colliding pair.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000)
        PARAMETER (MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sbe/nbe,nonbe,kbe(kszj,5),pbe(kszj,5),vbe(kszj,5)
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/sa26/ndiq(kszj),npt(kszj),ifcom(kszj),idi,idio
        common/sa28/nstr,nstra(kszj),nstrv(kszj),nstr0, &
         nstr1,nstr1a(kszj),nstr1v(kszj)
        common/trs/ntrs,nontrs,ktrs(kszj,5),ptrs(kszj,5),vtrs(kszj,5)


!       The hadronization model.
        i_had_model = INT( adj1(12) )
!       Moves parton ii (not equal to ik1 or ik2) to "trs".
        if( ( i_had_model == 0 .or. i_had_model == 3 ) .and. &
            ( i_proc == 6 .or. i_proc == 7 ) .and. &
            ( ii /= ik1 .or. ii /= ik2 ) )then
            ntrs = ntrs + 1
            ktrs(ntrs,1) = 3
            ktrs(ntrs,2) = k(ii,2)
            do j=1,4,1
                ptrs(ntrs,j) = P(ii,j)
                vtrs(ntrs,j) = V(ii,j)
            end do
            ptrs(ntrs,5) = P(ii,5)
            vtrs(ntrs,5) = 0D0
        end if

!       Updates diquark list.
        do j1=1,idi,1
            jj = ifcom(j1)
            kk = npt(j1)
            if(jj >= ii) ifcom(j1) = jj - 1
            if(kk >= ii) npt(j1)   = kk - 1
        end do
        do i2 = ii+1, N, 1
            i3 = i2 - 1
            ndiq(i3) = ndiq(i2)
        end do

!       Updates particle list "PYJETS".
        if(ii == N) N = N - 1
        if(ii < N)then
            do i1 = ii+1, N, 1
                i2 = i1 - 1
                do i3=1,5,1
                    K(i2,i3) = K(i1,i3)
                    P(i2,i3) = P(i1,i3)
                    V(i2,i3) = V(i1,i3)
                end do
                ishp(i2) = ishp(i1)
                taup(i2) = taup(i1)
            enddo
            N = N - 1
        end if
!       if(ik1 > ii) ik1 = ik1 - 1
!       if(ik2 > ii) ik2 = ik2 - 1

!       Updates particle list "sbe".
        if(ii == nbe) nbe = nbe - 1
        if(ii < nbe)then
            do i1 = ii+1, nbe, 1
                i2 = i1 - 1
                do i3=1,5,1
                    kbe(i2,i3) = kbe(i1,i3)
                    pbe(i2,i3) = pbe(i1,i3)
                    vbe(i2,i3) = vbe(i1,i3)
                end do
            end do
            nbe = nbe-1
        end if

!       Updates string list.
        do i1=1,nstr1,1
            nstraa = nstr1a(i1)
            if(nstraa >= ii) nstr1a(i1) = nstraa - 1
            nstrvv = nstr1v(i1)
            if(nstrvv >= ii) nstr1v(i1) = nstrvv - 1
        end do

!       Updates the values of lc(1-2,m) if which is  >=  ii.
        do m=1,icol,1
            lc1 = lc(1,m)
            if( lc1 >= ii ) lc(1,m) = lc1 - 1
            lc2 = lc(2,m)
            if( lc2 >= ii ) lc(2,m) = lc2 - 1
        end do

!       Removing colliding pairs with lc=0. This case happens when ii
!        is the 1-st one in PYJETS, especially when removing string
!        after inelastic collisions.
        jcol = 0
!       Loops over colliding pairs.
        do i=1,icol,1
            iic = lc(1,i)
            jjc = lc(2,i)
!           Throw away the pairs with lc=0.
            if(iic /= 0 .AND. jjc /= 0)then
                jcol = jcol + 1
                tc(1,jcol) = tc(1,i)
                tc(2,jcol) = tc(2,i)
                do k1=1,6,1
                    lc( k1, jcol ) = lc( k1, i )
                end do
                do k1=1,4,1
                    sig_tsmp_Cthetas_p( k1, jcol ) = sig_tsmp_Cthetas_p( k1, i )
                end do
            end if
        end do
        icol = jcol


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine strmov(istr1,istr1a,istr1v,ik1,ik2, i_proc )
!!      Removes all components of 'istr1'-th string out.
!       istr1: order number of string (in string list) to be moved
!       istr1a: line number (in parton list) of first component of above string
!       istr1v: line number (in parton list) of last component of above string
!       ik1 & ik2 are line # of colliding pair in parton list 'PYJETS'
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000)
        PARAMETER (MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sbe/nbe,nonbe,kbe(kszj,5),pbe(kszj,5),vbe(kszj,5)
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        common/sa26/ndiq(kszj),npt(kszj),ifcom(kszj),idi,idio
        common/sa28/nstr,nstra(kszj),nstrv(kszj),nstr0, &
         nstr1,nstr1a(kszj),nstr1v(kszj)
!       moves components of istr1-th string out of parton list


        do ii = istr1v, istr1a, -1
            call parmov( ii,ik1,ik2, i_proc )
        end do
!       moves string istr1-th out of string list
        if(istr1 == nstr1) nstr1 = nstr1 - 1
        if(istr1 < nstr1)then
            jj = istr1v - istr1a + 1
            do ii = istr1+1, nstr1, 1
!               if(ii >= istr1)then
!               nstr1v(ii-jj) = nstr1v(ii)
!               nstr1a(ii-jj) = nstr1a(ii)
                nstr1v(ii-1) = nstr1v(ii)
                nstr1a(ii-1) = nstr1a(ii)
!               end if
            end do
            nstr1 = nstr1 - 1
        end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine his_p( ic, jc, pi00, pj00, time_in, i_type )
!!      Propagates scattering particles to the collision time in the Lab frame.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        integer, intent(in) :: ic, jc, i_type
        real(kind=8), intent(in) :: pi00(4), pj00(4), time_in
        ! pi, pj: four-momenta before collisions
        ! i_type: event type. =1, elastic scattering; =2, inelastic; =3, shower.


        ! Zero formation time for scattered particles.

        ! Gives the formation time for the inelastic scattering (?).
        tau_p = 0D0
        if( i_type == 2 ) tau_p = 0D0
        time = time_in + tau_p
        do j=1,3,1
            V(ic,j) = V(ic,j) + pi00(j) / pi00(4) * ( time - V(ic,4) )
        end do
        V(ic,4)  = time
        taup(ic) = V(ic,4)

        ! Gives the formation time for the inelastic scattering (?).
        tau_p = 0D0
        if( i_type == 2 ) tau_p = 0D0
        time = time_in + tau_p
        do j=1,3,1
            V(jc,j) = V(jc,j) + pj00(j) / pj00(4) * ( time - V(jc,4) )
        end do
        V(jc,4)  = time
        taup(jc) = V(jc,4)


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine updpli_p( n00 )
!!      Updates collision time list for new partons when removing string
!!       and diqaruk after inelastic parton collisions.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/papr/t0,cspipiKK,dep,ddt,edipi,epin,ecsnn,ekn,ecspsn,ecspsm, &
         rnt,rnp,rao,rou0,vneu,vneum,ecsspn,ecsspm,ecsen
        common/scatt/ ic, jc, n0
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        integer, intent(in) :: n00
        integer :: l_coll(6)
        real(kind=8) :: t_coll, t_coll_min
        real(kind=8) :: sig_t_theta_p(4)
        logical :: succeed


!       Loops over new i (> n00) and old j partons (< n00).

        ! Time resolution.
        dddt = adj1(19)
        ! The collision pair counter (+1 avoids 0 icol).
        icol = icol + 1
        do i = n00+1, N, 1
            if( ishp(i) == 0 ) cycle
            ! Successful flag.
            succeed = .false.
            ! A number large enough to determine the minumum time.
            t_coll_min = 1D30
            ixc = 0
            jxc = 0
            loop_j: do j = 1, n00, 1
                if( ishp(j) == 0 ) cycle
                ! Avoids two particles colliding immediately after a collision.
                if( j == ic .OR. j == jc ) cycle
                call coij_p( i, j, i_fail, l_coll, t_coll, sig_t_theta_p )
                if( i_fail == 1 ) cycle
                ! Imposes the time resolution constraint.
                if( t_coll < 1D-10 ) cycle
                do j1 = 1, icol-1, 1
                    if( ABS( tc(1,j1) - t_coll ) < dddt ) cycle loop_j
                end do
                ! Chooses the smallest time for 'i'-cycle.
                if( t_coll < t_coll_min )then
                    lc(1,icol) = l_coll(1)
                    lc(2,icol) = l_coll(2)
                    lc(3,icol) = l_coll(3)
                    lc(4,icol) = l_coll(4)
                    lc(5,icol) = l_coll(5)
                    lc(6,icol) = l_coll(6)
                    tc(1,icol) = t_coll
                    tc(2,icol) = t_coll
                    sig_tsmp_Cthetas_p(1,icol) = sig_t_theta_p(1)
                    sig_tsmp_Cthetas_p(2,icol) = sig_t_theta_p(2)
                    sig_tsmp_Cthetas_p(3,icol) = sig_t_theta_p(3)
                    sig_tsmp_Cthetas_p(4,icol) = sig_t_theta_p(4)
                    t_coll_min = t_coll
                    ixc = l_coll(1)
                    jxc = l_coll(2)
                end if
                succeed = .true.
            end do loop_j
            ! Keeps the one with the smallest time from pairs including i or j.
            if( succeed )then
                n_jump_out = 0
                j1 = 1
                do while(.true.)
                    if( j1 > icol-1 ) exit
                    iic = lc(1,j1)
                    jjc = lc(2,j1)
                    if(       ixc /= iic .AND. ixc /= jjc &
                        .AND. jxc /= iic .AND. jxc /= jjc )then
                        j1  = j1 + 1
                        cycle
                    end if
                    ttc = tc(1,j1)
                    ! Thows away the pair with larger time.
                    if( ttc > t_coll_min )then
                        k_begin = j1
                        n_jump_out = n_jump_out + 1
                    else
                        k_begin = icol
                        n_jump_out = 2
                    end if
                    icol = icol - 1
                    do k1 = k_begin, icol, 1
                        lc( 1, k1 ) = lc( 1, k1+1 )
                        lc( 2, k1 ) = lc( 2, k1+1 )
                        lc( 3, k1 ) = lc( 3, k1+1 )
                        lc( 4, k1 ) = lc( 4, k1+1 )
                        lc( 5, k1 ) = lc( 5, k1+1 )
                        lc( 6, k1 ) = lc( 6, k1+1 )
                        tc( 1, k1 ) = tc( 1, k1+1 )
                        tc( 2, k1 ) = tc( 2, k1+1 )
                        sig_tsmp_Cthetas_p(1,k1) = sig_tsmp_Cthetas_p( 1, k1+1 )
                        sig_tsmp_Cthetas_p(2,k1) = sig_tsmp_Cthetas_p( 2, k1+1 )
                        sig_tsmp_Cthetas_p(3,k1) = sig_tsmp_Cthetas_p( 3, k1+1 )
                        sig_tsmp_Cthetas_p(4,k1) = sig_tsmp_Cthetas_p( 4, k1+1 )
                    end do
                    if( n_jump_out == 2 ) exit
                end do
                ! For the next time finding.
                icol = icol + 1
            end if
        end do
        ! Deducts 1 because the counter begins with +1 not 0.
        icol = icol - 1
        n0 = N


        return
        end



!ccccccccccccccccccccccccccccccccccccc end ccccccccccccccccccccccccccccccccccccc